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Abstract 

Computational  fluid  dynamics  (CFD)  methods  have  been  coupled  with  structural 
solvers  to  provide  accurate  predictions  of  limit  cycle  oscillations  (LCO).  There  is,  however, 
a  growing  interest  in  understanding  how  uncertainties  in  flight  conditions  and  structural 
parameters  affect  the  character  of  an  LCO  response,  leading  to  failure  of  an  aeroelas- 
tic  system.  Uncertainty  quantification  of  a  stochastic  system  (parametric  uncertainty) 
with  stochastic  inputs  (initial  condition  uncertainty)  has  traditionally  been  analyzed 
with  Monte  Carlo  simulations  (MCS).  Probability  density  functions  (PDF)  of  the  LCO 
response  are  obtained  from  the  MCS  to  estimate  the  probability  of  failure.  A  CFD  so¬ 
lution,  however,  can  take  days  to  weeks  to  obtain  a  single  response,  making  the  MCS 
method  intractable  for  large  problems.  A  candidate  approach  to  efficiently  estimate  the 
PDF  of  an  LCO  response  is  the  stochastic  projection  method.  The  classical  stochastic 
projection  method  is  a  polynomial  chaos  expansion  (PCE).  The  PCE  approximates  the  re¬ 
sponse  in  the  stochastic  domain  through  a  Fourier  type  expansion  of  the  Wiener-Hermite 
polynomials.  An  LCO  response  can  be  characterized  as  a  subcritical  or  supercritical  bi¬ 
furcation,  and  bifurcations  are  shown  to  be  discontinuities  in  the  stochastic  domain.  The 
PCE  method,  then,  would  be  too  inefficient  for  estimating  the  LCO  response  surface. 
The  objective  of  this  research  is  to  extend  the  stochastic  projection  method  to  include 
the  construction  of  B-spline  surfaces  in  the  stochastic  domain.  The  multivariate  B-splinc 
problem  is  solved  to  estimate  the  LCO  response  surface.  An  MCS  is  performed  on  this 
response  surface  to  estimate  the  PDF  of  the  LCO  response.  The  probability  of  failure 
is  then  computed  from  the  PDF.  The  stochastic  projection  method  via  B-splines  is  ap¬ 
plied  to  the  problem  of  estimating  the  PDF  of  a  subcritical  LCO  response  of  a  nonlinear 
airfoil  in  inviscid  transonic  flow.  The  stochastic  algorithm  provides  a  conservative  esti¬ 
mate  of  the  probability  of  failure  of  this  aeroelastic  system  two  orders  of  magnitude  more 
efficiently  than  performing  an  MCS  on  the  governing  equations. 
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QUANTIFYING  INITIAL  CONDITION  AND  PARAMETRIC 
UNCERTAINTIES  IN  A  NONLINEAR  AEROELASTIC  SYSTEM 
WITH  AN  EFFICIENT  STOCHASTIC  ALGORITHM 


I.  Introduction 

Fluid-structure  interaction  can  result  in  the  loss  of  dynamic  stability  to  a  time  periodic 
instability  that  grows  unbounded,  i.e.,  flutter  [1],  When  nonlinear  aerodynamics  (e.g., 
transient  shocks  or  boundary  layer  separation)  or  nonlinear  structural  parameters  are 
present  to  counter  the  growth  of  the  unstable  mode,  the  dynamic  response  may  stabilize 
to  a  limit  cycle  oscillation  (LCO)  [2],  The  study  of  LCO  in  aeroelastic  systems  is  an  active 
area  of  research  [3,  4,  5,  6].  Flight  data  from  fighter  aircraft  with  external  stores  indicate 
that  the  LCO  response  is  characterized  by  antisymmetric  motion  of  the  wing  and  a  lateral 
motion  of  the  fuselage  and,  consequently,  the  aircrew.  LCO  can  lead  to  fatigue  failure  of 
the  wing  structure.  The  lateral  motion  LCO  induces  can  degrade  mission  effectiveness  by 
not  allowing  the  aircrew  to  track  a  target  or  lead  the  aircrew  to  believe  that  the  aircraft 
has  begun  to  flutter,  causing  a  mission  abort.  In  severe  cases,  LCO  represents  a  hazard 
to  flight  safety  [2], 

Computational  fluid  dynamics  (CFD)  methods  have  been  coupled  with  structural 
solvers  to  provide  accurate  predictions  of  LCOs  [7].  There  is,  however,  a  growing  interest 
in  understanding  how  uncertainties  in  flight  conditions  and  structural  parameters  affect 
the  character  of  an  LCO  response,  leading  to  failure  of  an  aeroelastic  system  [8].  Un- 
certainty  quantification  of  a  stochastic  system  (parametric  uncertainty)  with  stochastic 
inputs  (initial  condition  uncertainty)  has  traditionally  been  analyzed  with  Monte  Carlo 
simulations  (MCS)  [9].  In  an  MCS,  N  random  realizations  of  the  uncertainty  parameters 
are  input  into  the  system  of  equations  to  obtain  N  realizations  of  the  response.  The 
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probability  density  function  (PDF)  of  the  response  is  then  approximated  from  the  N  re¬ 
alizations.  The  MCS  method  only  requires  that  the  distribution  of  the  input  parameters 
be  known.  A  CFD  solution,  however,  can  take  days  to  weeks  to  obtain  a  single  response, 
making  the  MCS  method  intractable  for  large  problems. 

Other  techniques  that  have  been  used  to  solve  stochastic  systems  with  stochastic 
inputs  are  based  upon  Galerkin’s  method  [9].  These  include  reduced  order  modeling 
(ROM)  [10,  11,  12]  and  polynomial  chaos  expansions  (PCE)  [13,  14,  15,  9,  16,  17,  18, 
19,  20,  21].  The  reduced  order  model  represents  a  full  order  solution  with  an  optimal 
basis,  in  the  mean  square  sense,  through  a  Karhunen-Loeve  expansion.  Monte  Carlo 
simulations  can  efficiently  be  run  on  the  ROM,  but  only  for  small  variations  in  the  inputs. 
Thus,  ROM  is  useful  for  examining  responses  near  the  mean  value  but  cannot  predict 
responses  with  large  variations  in  the  input  [10].  The  PCE  method,  on  the  other  hand, 
allows  large  variations  in  the  input  parameters  and  has  been  shown,  for  certain  classes 
of  problems,  to  be  far  more  efficient  than  an  MCS  [15].  The  classical  PCE  method,  also 
known  as  the  Wiener- Hermite  expansion  [22],  falls  into  a  general  category  of  stochastic 
projection  methods.  Of  the  two  methods,  ROM  or  stochastic  projection,  the  stochastic 
projection  method  was  chosen,  clue  to  its  greater  efficiency  and  robustness  in  handling 
large  variations,  to  quantify  the  effects  of  uncertainty  and  estimate  the  probability  of 
failure  for  a  nonlinear  aeroclastic  system. 

1.1  Overview  of  Stochastic  Projection  Methods 

Stochastic  projection  methods,  unlike  the  Karhunen-Loeve  expansions  employed  by 
ROM,  provide  a  way  to  estimate  the  expansion  of  a  response  without  a  priori  knowledge 
of  the  form  of  the  response  in  the  stochastic  domain.  A  random  basis,  denoted  \&j(£), 
orthogonal  with  respect  to  the  distribution  of  the  input  uncertainty,  is  typically  selected. 
The  deterministic  coefficients  of  the  expansion,  cq(f),  are  either  computed  as  part  of  the 
time  integration  (intrusive  approach)  or  estimated  from  a  limited  MCS  (non-intrusive 
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approach).  The  expansion  of  some  response  a(£,  £)  then  takes  the  form  [14,  15] 


p 

«(*,£)  (i-1) 

i= 0 

where  the  upper  limit  P  is  based  on  the  order  of  the  expansion  and  the  number  of  input 
uncertainties,  i.e.  the  dimension,  d,  of  The  vector  £  contains  the  random  variables 
that  characterize  the  uncertainty  distribution  of  the  input  parameters.  For  a  Gaussian 
distribution,  the  uncertainty  of  some  input  parameter,  (3a,  is  characterized  by  its  mean 
value,  (3a,  and  its  standard  deviation,  /3a,  as 

Pa  —  Pa  +  £iA*  (1-2) 


where  is  a  zero  mean,  unit  variance  Gaussian  random  variable. 

With  the  intrusive  approach,  the  expansion  in  Eq.  (1.1)  is  substituted  directly  into 
the  governing  equations  of  motion.  For  example,  if  an  equation  of  motion  took  the  form 

C[a(t)}  +  paa3  =  0  (1.3) 


where  £  is  a  linear  differential  operator,  then  the  stochastic  projection  method  would, 
with  a  Galerkin  approach  (see,  for  example,  Le  Maitre  et  al.  [18]),  produce  the  following 
system  of  equations  to  be  solved 


rtz,  u\-\  ,  a  \  \  \  -  TI'/T./T/,.  vlG)  ,  Q 

/  J/  J/  jOLjOLjOLk  lTt  Tt  v 


i= 0  j= 0  k= 0 


(*n,  *n) 


\  ^  \  ^  \  ^  A  ~  V|/n)  _  n 

/  J  /  J  /  J  T  ~w~  \  0 


2=0  j= 0  k= 0 


(*n,  *n) 


(1.4) 


where  the  free  index,  n,  is  in  the  range  0  <  n  <  P.  The  expected  value  operator,  (  ),  is 
defined  as 


r»oo  /»oo 


</(*),*»(*)>  = 


d^ldG-'-d^ 


(1.5) 


'  — oo  J  — oo  «/  —  oo 
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where  /(£)  is  any  function  of  £  and  w(£)  is  the  Guassian  PDF 


“>(£) 


(1.6) 


The  condition  of  orthogonality  with  respect  to  the  Guassian  PDF  substantially  reduces 
the  number  of  non-zero  expected  values  for  a  linear  problem.  As  nonlinearity  increases, 
the  number  of  non-zero  terms  grows  rapidly  due  to  the  coupling  evident  in  Eq.  (1.4). 
Computational  efficiency  is  rapidly  lost  when  high  order  expansions  of  highly  nonlinear 
equations  are  required  [22,  23,  24,  25,  26]. 

To  gain  back  some  computational  efficiency,  researchers  have  employed  the  non- 
intrusive  approach.  This  is  a  sampling  approach  and  considered  non-intrusive  since 
parametric  uncertainty  is  not  included  analytically  in  the  equations  of  motion.  A  limited 
MGS  is  performed  and  the  responses  obtained  are  used  to  estimate  the  coefficients  of  the 
expansion.  Due  to  the  random  selection  of  samples,  the  expected  value  operator  is  again 
employed.  The  coefficients  are  estimated  from  [8] 


(Mh) 

('I'n,  *n>  ' 


(1.7) 


This  non-intrusive  approach  substantially  reduces  the  number  of  expected  values  that 
need  to  be  pre-computed  and  stored  and  is  considered  more  appropriate  for  highly  non¬ 
linear  problems. 

Thus  far  the  discussion  has  centered  around  the  mathematics  of  the  stochastic 
projection  method.  The  stochastic  projection  also  has  a  physical  interpretation,  i.e.  it 
can  be  viewed  as  approximating  the  response  surface  across  the  entire  stochastic  domain 
[9,  15,  25].  The  stochastic  projection  method,  though,  differs  from  traditional  response 
surface  methodologies  (RSM)  [27]  in  that  the  expansion  is  global  whereas  RSM  is  con¬ 
cerned  with  local  expansions  about  some  point  of  interest.  Thus,  the  stochastic  projection 
method  is  an  interpolation  algorithm.  This  viewpoint  explains  many  of  the  issues  that, 
in  particular,  caused  the  PCE  approach  to  be  mostly  abandoned  by  the  1970’s  [15].  The 
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original  PCE  method  assumed  a  Gaussian  uncertainty  in  the  parameters  and  the  orthog¬ 
onal  basis  functions  chosen  for  the  expansion  were  the  Hermite  polynomials  [13,  14,  15]. 
Chorin  [22]  discussed  the  lack  of  sufficient  frequency  content  in  the  Hermite  polynomials 
to  approximate  discontinuities  to  stochastic  solutions  of  the  inviscid  Burger’s  equation. 
Crow  and  Canavan  [23]  demonstrated  that,  while  only  low  order  Hermite  expansions 
are  computationally  practical,  high  order  Hermite  polynomials  are  required  to  allow  an 
energy  cascade  for  approximating  a  turbulent  flow. 

Despite  these  findings,  recent  literature  on  the  generalized  PCE  through  an  Askey 
scheme  [19],  which  matches  distribution  functions  to  the  appropriate  orthogonal  poly¬ 
nomials,  still  cites  the  optimality  of  polynomial  expansions  without  regard  to  the  issue 
of  frequency  content.  Millman  et  al.  [24,  25],  however,  showed  that  by  changing  to  a 
Fourier  basis,  termed  a  Fourier  chaos  expansion  (FCE),  stochastic  projections  with  dis¬ 
continuities  could  be  better  approximated.  A  comparison  of  the  stochastic  projections 
from  the  PCE  and  FCE  approaches  is  shown  in  Figure  1.1,  where  for  this  example  the 
random  variable  £2  represents  the  uncertainty  in  the  initial  pitch  angle  of  a  nonlinear 
airfoil  in  the  presence  of  modeled  aerodynamics  [25].  The  ability  of  the  FCE  to  better 
approximate  the  discontinuity  than  the  PCE  method  leads  to  a  better  approximation  of 
the  PDF  of  the  LCO  response,  as  shown  in  Figure  1.2.  These  results  were  based  on  the 
intrusive  method  and,  while  a  substantial  improvement  over  the  PCE  method,  the  FCE 
method  was  still  too  computationally  inefficient  to  be  useful  in  a  design  environment  due 
to  multiple  summations  such  as  the  ones  in  Eq.  (1.4). 

Numerical  techniques  have  been  advanced  to  overcome  the  computational  expense 
of  computing  the  multiple  summations.  For  example,  Mathelin  and  Hussaini  [28]  employ 
a  stochastic  collocation  scheme  in  the  solution  of  the  Reimann  problem  that  reduces  the 
computational  effort  from  order  P3,  0(P3),  to  O(P).  However,  they  note  that  the  upper 
limit  P  has  to  be  large  (P  =  200)  to  accurately  account  for  the  moving  discontinuity 
[28]- 
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-  Monte  Carlo 

- 8th  order  PCE 


-  Monte  Carlo 

- 8,h  order  FCE 


(a)  MCS  vs.  8th  order  PCE 


(b)  MCS  vs.  8th  order  FCE 


Figure  1.1  a(£,  £)  at  t  =  2000  vs.  £2  [25] 


Figure  1.2  PDFs  of  a  supercritical  response,  MCS  vs.  8th  order  PCE  and  FCE  [25] 


It  will  be  shown  that  bifurcations,  such  as  the  LCO  response  of  an  airfoil,  are 
discontinuities  in  the  stochastic  domain.  Thus,  the  numerical  issues  faced  by  Chorin  [22] , 
Crow  and  Canavan  [23],  and  Mathelin  and  Hussaini  [28]  must  be  dealt  with  in  an  efficient 
manner.  In  interpolation  theory,  discontinuities  are  treated  as  piecewise  polynomials  and 
researchers  commonly  employ  B-splines  to  approximate  these  piecewise  polynomials  [29] . 
Constructing  B-spline  surfaces  is  also  a  collocation  approach  since  the  nodes  of  the  surface 
are  fitted  exactly  [30].  The  nodes  of  the  response  (B-spline)  surface  can  be  computed 
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from  a  nonintrusive  approach  and,  since  the  B-spline  method  is  non-Galerkin,  expected 
values  need  not  be  computed  or  stored  [31]. 

1.2  Scope  of  Research 

The  objective  of  this  research  is  to  extend  the  stochastic  projection  method  to 
include  the  construction  of  B-spline  surfaces  in  the  stochastic  domain.  These  surfaces 
should  adequately  resolve  discontinuities  that  arise  from  the  bifurcations  (LCO  response) 
of  a  nonlinear  airfoil  in  an  inviscid  transonic  flow.  An  efficient  MCS  over  the  response 
surface  will  allow  a  rapid  estimation  of  the  PDF  of  an  LCO  response.  The  probability 
of  failure  can  then  be  computed  from  the  PDF.  The  scope  of  this  research  is  defined  by 
the  following  thesis  statement,  approach,  and  assumptions: 

1.2.1  Thesis  Statement.  The  probability  density  function  of  a  limit  cycle  os¬ 
cillation  response,  and  hence  the  probability  of  failure,  can  be  accurately  and  efficiently 
estimated  with  a  small  sample  size  obtained  from  the  solution  of  the  nonlinear  time  de¬ 
pendent  aeroelastic  system  of  equations. 

1.2.2  Research  Approach  and  Assumptions.  The  implementation  and  effec¬ 
tiveness  of  the  stochastic  projection  method  via  B-splines  will  be  demonstrated  for  the 
CFD  application  of  a  two-dimensional  pitching  and  plunging  airfoil  in  inviscid,  transonic 
compressible  flow.  The  structural  components  will  be  nonlinear  in  pitch.  Since  an  MCS 
is  impractical  for  the  CFD  application,  the  stochastic  projection  method  via  B-splines 
will  first  be  validated  against  MCS  results  of  a  model  problem:  a  structurally  nonlinear 
airfoil  in  a  linearly  modeled,  low  speed  incompressible  flow.  Convergence  of  the  PDFs 
for  the  CFD  application  will  then  be  presented. 

The  uncertainty  distributions  in  this  research  effort  are  limited  to  Gaussian  distri¬ 
butions,  as  given  in  Eq.  (1.6).  The  choice  of  Gaussian  distributions  was  made  early  in  the 
research  phase  to  allow  a  direct  comparison  with  the  classical  PCE  formulation.  As  will 
be  shown,  the  choice  in  distributions  is  not  critical  to  the  development  of  this  stochastic 
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algorithm.  The  number  of  uncertainties  was  limited  to  one  initial  condition  uncertainty 
and  one  parametric  uncertainty. 

1.3  Document  Organization 

The  remainder  of  this  document  is  organized  as  follows: 

Chapter  II:  Demonstrates  estimating  the  probability  of  failure  of  a  nonlinear  airfoil 
in  flow  modeled  with  linear  aerodynamics  through  a  traditional  Monte  Carlo  analysis. 

Chapter  III:  Develops  the  theory  of  the  stochastic  projection  method  via  B-splines. 
The  stochastic  algorithm  is  implemented  on  the  nonlinear  airfoil  with  linear  aerodynamics 
to  estimate  the  probability  of  failure. 

Chapter  IV:  Describes  the  development  of  an  inviscid  aeroclastic  code  (EULER- 
AE)  that  computes  the  flow  properties  over  a  pitching  and  plunging  airfoil  with  the  Euler 
equations. 

Chapter  V:  Implements  the  stochastic  algorithm  on  the  nonlinear  airfoil  in  inviscid 
transonic  compressible  flow  to  estimate  the  probability  of  failure. 

Chapter  VI:  Summarizes  the  results  and  conclusions  of  the  current  research.  Dis¬ 
cusses  topics  of  future  research  suggested  by  the  currrent  research. 
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II.  Nonlinear  Pitch  and  Plunge  Airfoil  with  Linear  Aerodynamics 

2. 1  Introduction 

Before  choosing  an  expansion  to  represent  a  response  in  the  stochastic  domain,  it  is 
useful  to  examine  the  behavior  of  a  set  of  nonlinear  equations  with  parametric  and  initial 
condition  uncertainties.  An  appropriate  model  problem  to  begin  this  analysis  should  be 
easily  integrated  in  time  such  that  an  MCS  is  practical.  To  this  end,  a  nonlinear  pitch 
and  plunge  airfoil  with  linear  aerodynamics  is  chosen.  A  stability  analysis  of  this  model 
problem  was  accomplished  by  Lee  et  al  [3].  The  goal  of  this  chapter  is  not  to  replicate 
the  stability  analysis  of  a  system  of  ordinary  differential  equations  with  the  methods  one 
finds  in,  for  example,  Seydel  [32],  but  rather,  to  perform  a  numerical  analysis  of  the 
behavior  of  the  bifurcations  of  a  system  in  the  presence  of  uncertainties.  Results  from 
Lee  et  al.  [3]  and  Beran  and  Pettit  [8]  are  a  basis  for  the  validity  of  the  numerical  analysis 
performed  here. 

The  equations  of  motion  for  a  nonlinear  pitch  and  plunge  airfoil  are  first  derived. 
The  aerodynamics  needed  to  compute  the  lift  and  moment  coefficients  are  obtained  from 
Fung’s  model  [1]  for  incompressible  flow  based  on  the  circulation  theory  of  lift.  A  set 
of  ordinary  differential  equations  are  then  derived,  following  Lee  et  al  [3].  PDFs  of  the 
LCO  response  in  pitch,  aLco,  are  developed  from  an  MCS  based  on  parametric  and 
initial  condition  uncertainty.  These  PDFs  are  also  used  to  generate  the  probability  of 
failure  based  on  encountering  an  LCO  response.  This  analysis  will  serve  as  basis  for 
the  development  in  the  next  chapter  of  a  stochastic  algorithm  that  will  approximate  the 
PDFs  obtained  here  without  performing  an  MCS  on  the  governing  equations. 

2.2  Equations  of  Motion 

Various  formulations  for  the  pitching  and  plunging  airfoil  have  been  developed 
[1,  3,  33,  34],  The  notation  shown  in  Figure  2.1  is  used  for  the  derivation  of  the  equations 
of  motion  for  this  model  problem  [1,  3].  The  plunge,  h,  is  measured  positive  in  the 
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Figure  2.1  Schematic  of  a  pitch  and  plunge  airfoil 


downward  direction.  The  elastic  axis  (or  hinge  point)  is  located  the  distance  cihb  from 
the  midchord,  b,  and  is  measured  positive  aft.  The  pitch  angle,  a ,  is  positive  nose  up 
about  the  elastic  axis.  The  center  of  mass  is  located  the  distance  xab  from  the  elastic 
axis,  also  measured  positive  aft.  Summing  the  forces  and  moments,  and  assuming  small 
pitch  angles,  results  in 


mh  +  Saa  +  Khh  —  —L  (2.1a) 

Sah  T  Ia&  T  T  /3Qa3  T  7aQ^)  =  -Tfea  (2.1b) 

where  L  is  the  lift,  Mea  is  the  moment  about  the  elastic  axis,  Sa  is  the  first  moment  of 
inertia,  I0  is  the  second  moment  of  inertia,  Kh  is  the  plunge  spring  constant,  and  Ka 
is  the  pitch  spring  constant.  The  constants  (3a  and  7 a  control  the  amount  of  structural 
nonlinearity.  Note  that  [1] 


=  mxab 

(2.2a) 

la 

=  mrab 

(2.2b) 

Kh 

=  mul 

(2.2c) 

Ka 

= 

(2. 2d) 
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where  m  is  the  mass  of  the  airfoil,  rab  is  the  radius  of  gyration,  and  u>h  and  uia  are  the 
uncoupled  natural  frequencies  of  plunge  and  pitch,  respectively.  Also 


L  =  TiPooV^cQ 


2 

1 


Mea  =  -PooV^c  C„ 


(2.3a) 

(2.3b) 


where  c  is  the  chord  length,  p is  the  freestream  density,  VA,  is  the  freestream  velocity, 
and  Ci  and  Cm  are  the  lift  and  moment  coefficients,  respectively.  By  defining  the  following 
nondimensional  parameters 


(2.4a) 

(2.4b) 

(2.4c) 

(2.4d) 


h 

y  = 

b 

Uh 

uy  = 

— 

Vr  = 

Voo 

ujab 

m 
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Poob2ir' 


where  y  is  the  nondimensional  plunge,  ay  is  the  frequency  ratio,  Vr  is  the  reduced  velocity, 
and  /is  is  the  mass  ratio,  and  differentiating  with  respect  to  the  nondimensional  time 


tv« 


T  = 


(2.5) 


Eqs.  (2.1a)  and  (2.1b)  can  be  rewritten  as 


ay 


y  +  xaa  +  (  —  I  y  = - Ct 

.  Vr  )  nps 

x  1  2 

—y  +  ol  +  7yy(«  +  (3aot 3  +  ba®5)  = - 

ri  V 


(2.6a) 

(2.6b) 


Both  the  lift  and  moment  coefficients  are  a  function  of  time.  Lee  et  al.  [3]  provide 
expressions  for  incompressible  flow  over  a  pitching  and  plunging  airfoil  based  on  a  model 
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given  by  Fung  [1] 


C/(t)  =  7r  (y  —  aha  +  a)  +  2n 


+  27 r  /  0(r  —  a) 


a(0)  +7/(0)  +  (  -  -  ah  )  d(0) 


dcr 


0(r) 


d(a)  +  y(a)  +  (  -  -  ah  )  d(a) 


(2.7a) 


Cm(r)  = 


=  7T  |  -  +  a* 


a(0)  +  3/(0)  +  (  -  -  ah  )  d(0) 


7T  I  -  +  ah  )  /  4>{t  —  a) 


0(r) 


d(c)  +  j/(c)  +  l  -  -  a,h  J  afaj 


da 


(2.7b) 


+  )7ra*  (i/  -  ah&)  a,M  ®  - 


where  the  Wagner  function,  0(r),  is  given  by  [1] 


4>{t)  =  1  -  0i  e  £1T  -  -02  e 


— £2  T 


(2.8) 


Fung  [1]  provides  several  values  for  the  constants  in  the  Wagner  function.  The  values  of 
=  0.165,  02  =  0.335,  £\  =  0.0455,  and  e2  =  0.3  are  obtained  from  Jones  [35]. 


2.3  Numerical  Integration 

Equations  (2.6a),  (2.6b),  (2.7a)  and  (2.7b)  represent  a  complete  description  of  a 
pitching  and  plunging  airfoil  with  incompressible  aerodynamics.  Since  stability  analysis 
techniques  require  a  system  of  ordinary  differential  equations,  Lee  et  al.  [3]  introduced 
four  new  variables 


wi  (r) 


u>2(t) 


W3{t) 


e~£1 (T" ff)a(a)  da 

g— £2(T  — <T  )a(a) 
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(2.9a) 

(2.9b) 

(2.9c) 
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w4(t)  —  f  e  e2<T  aPi(a)  cl  a.  (2.9d) 

Jo 

These  variables  allow  Eqs.  (2.6a)  and  (2.6b)  to  be  rewritten  as 


c0y  +  cib  +  c2y  +  c3a  +  cAy  +  c5a  +  c6w1  +  c7w2  +  c$w3  +  c9wA  =  f(r) 
d0y  T  d\di.  -t-  d2oi  T  d3y  T  d^cx  T  d3c\J  T  d^oP  T  d7y 

+d8w  1  +  d9w2  +  dww3  +  diiW4  =  g(r) 


(2.10a) 

(2.10b) 


where  the  constants  in  Eq.  (2.10a)  are  given  by 
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and  the  constants  in  Eq.  (2.10b)  are  given  by 
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dg 
d  10 
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The  forcing  terms,  f(r)  and  g{r)  are  given  by 


I(t)  = 


2 


-  ah  a(0)  +  2/(0) 


9(T)  = 


1  +  2  ah. 

K 


(tpiEie  slT  +  4r2£2e  S2T)  (2.11a) 

/(r).  (2.11b) 


With  this  formulation,  the  equations  of  motion  can  be  separated  into  eight  first  order, 
ordinary  differential  equations.  By  letting  xA  =  a,  x2  =  d,  x3  =  y,  x A  =  y,  x5  —  Wi, 
Xq  =  w2,  x7  =  w3,  and  Xg  =  W4,  the  following  equations  result 


■i\ 

=  x2 

(2.12a) 

X2 

coH  —  d0P 

cid0  —  CqcIi 

(2.12b) 

=  X4 

(2.12c) 

XA 

-e  8 

1  1 

0  0 

(2.12d) 

x5 

=  XA~  £AX5 

(2.12e) 

xe 

=  Xi  -  £2X6 

(2.12f) 

x7 

=  x3-  £1X7 

(2.12g) 

Xg 

=  X3  ~  £2X8 

(2.12h) 

where 


P  =  c2xA  +  c3x2  +  cAx3  +  c5Xi  +  c6x 5  +  c7x6  +  c8x7  +  c9x8  -  f(r)  (2.13a) 
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and 


H  =  d2X  2  +  d3X4  +  di^xi  +  d$x\  +  d§x\  +  ^7X3 


+  d%x  5  +  dgXQ  +  dio^7  +  dnx$  —  d{T)- 


(2.13b) 


With  the  aeroelastic  equations  formulated  as  a  system  of  first  order,  ordinary  dif¬ 
ferential  equations,  a  variety  of  methods  are  available  for  numerical  integration.  The 
method  chosen  here  is  a  simple  forward  Euler  integration.  The  time  step  should  be  1/256 
of  the  shorter  period  of  the  two  coupled  modes  of  oscillation  for  an  accurate  time  inte¬ 
gration  [3].  Lee  et  al.  [3]  used  a  time  step  of  At  =  0.1  to  satisfy  this  requirement.  The 
same  time  step  was  used  in  this  study. 

2-4  Bifurcation  Behavior 

Bifurcations  are  a  change  in  the  stability  of  a  system  that  results  in  different  dy¬ 
namic  responses  [32],  A  linear  aeroelastic  system,  at  some  critical  reduced  velocity,  will 
undergo  a  change  from  a  stable  stationary  (no  motion)  solution  to  divergent  flutter  (un¬ 
bounded  motion).  This  critical  point  is  termed  by  aeroelasticians  as  the  flutter  point 
[1],  Nonlinear  aeroelastic  systems  can  stabilize  beyond  the  flutter  point  to  a  periodic 
dynamic  response  [3,  2,  4],  Mathematicians  term  this  change  from  a  stationary  response 
to  a  periodic  response  the  Hopf  bifurcation  point  [32],  Note  that  the  flutter  point  and 
the  Hopf  bifurcation  point  are  coincident  for  an  aeroelastic  system.  The  periodic  motion 
beyond  the  flutter  point  has  been  termed  limited  amplitude  flutter  or  limit  cycle  oscilla¬ 
tion  (LCO).  While  the  latter  term  is  more  a  mathematical  term,  its  use  is  prevalent  in 
the  aeroelastic  literature  [2],  In  this  study,  the  terms  flutter  point  and  LCO  are  used. 

Two  different  types  of  dynamic  responses  are  examined  in  this  study.  The  first  type 
of  response  was  described  above  and  is  termed  a  supercritical  response.  A  supercritical 
response  is  a  stable  periodic  solution  for  reduced  velocity  beyond  the  flutter  point.  As 
the  reduced  velocity  decreases  through  the  flutter  point,  the  periodic  solution  reduces  to 
the  stationary  solution.  A  supercritical  response  was  obtained  by  Beran  and  Pettit  [8] 
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and  Millman  et  al.  [31]  from  the  aeroclastic  system  described  in  the  previous  sections 
with  the  following  structural  parameters: 

jis  =  100  ah  =  —0.5 

u or  =  0.2  Pot  =  3 

la  =  20 

The  value  for  7Q  was  selected  to  to  stabilize  a  subcritical  response,  which  will  be  discussed 
shortly.  For  this  supercritical  case,  both  Beran  and  Pettit  [8]  and  Millman  et  al.  [31] 
located  the  flutter  point  at  Vr  ~  6.28.  The  time  histories  of  pitch  at  Vr  =  6.1,  6.2,  6.3 
and  6.4  with  an  initial  pitch  of  a(0)  =  1°  are  shown  in  Figure  2.2.  The  remaining  initial 
conditions  are  given  by  d(0)  =  y( 0)  =  y(0)  =  0.  For  reduced  velocities  below  the  flutter 
point,  Figs.  2.2(a)  and  (b),  the  initial  pitch  displacement  damps  to  no  motion.  At  reduced 
velocities  above  the  flutter  point,  Figs.  2.2(c)  and  (d),  the  initial  pitch  displacement 
causes  the  airfoil  to  exhibit  an  LCO.  Also  note  that  a  longer  time  integration  is  required 
to  obtain  the  final  motion  of  the  airfoil  the  closer  the  reduced  velocity  is  to  the  flutter 
point.  By  plotting  the  LCO  amplitude  at  various  reduced  velocities,  the  bifurcation 
diagram  in  Figure  2.3  results.  Time  integrations  near  the  flutter  point  were  carried  out 
to  r  =  20,  000  before  the  stable  response  (stationary  or  LCO)  was  achieved. 

Lee  et  al.  [3]  demonstrated  that  a  subcritical  response,  i.e.,  a  dynamically  unstable 
response,  results  when  a  negative  value  for  the  cubic  spring  constant,  /3a,  and  a  zero 
quintic  restoring  force  (ya  =  0)  are  applied  [3].  At  reduced  velocities  above  the  flutter 
point,  the  airfoil  will  exhibit  divergent  flutter  for  any  initial  pitch  except  a(0)  =  0,  the 
dynamically  unstable  stationary  branch.  With  a  large  initial  pitch,  divergent  flutter  can 
also  be  exhibited  at  reduced  velocities  below  the  flutter  point.  A  bifurcation  diagram  of 
a  subcritical  response  is  shown  in  Figure  2.4.  The  source  of  this  instability  can  be  seen 
in  the  pitch  spring  term 

Act(o:  +  Pao?  +  7aci5) 


xa  =  0.25 
ra  =  -0.5 
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a  (degrees)  a  (degrees) 


(a)  Vr  =  6.1  (b)  Vr  =  6.2 


Figure  2.2  Time  histories  of  pitch  with  a(0)  =  1.0° 


Figure  2.3  Bifurcation  diagram  of  a  supercritical  response 
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Figure  2.4  Bifurcation  diagram  of  a  subcritical  response 

in  Eq.  (2.1b).  With  the  absence  of  the  quintic  spring  constant  7„,  a  negative  value  for 
cubic  spring  constant  f3a  can  be  viewed  as  loosening,  or  weakening,  the  pitch  spring.  For 
sufficiently  large  values  of  pitch,  the  negative  cubic  spring  constant  dominates  the  pitch 
spring  resulting  in  a  loss  of  stability. 

In  order  to  regain  stability  in  the  structure,  Millman  et  al.  [25]  added  the  quintic 
restoring  force  in  Eq.  (2.1b).  At  large  amplitudes,  the  quintic  spring  constant  7a  domi¬ 
nates  the  pitch  spring  and  the  structure  exhibits  a  large  amplitude  LCO.  Following  Lee 
et  al.  [3],  the  cubic  spring  constant  was  changed  from  f3a  =  3  to  f3a  =  —3.  Time  histories 
of  pitch  are  shown  in  Figure  2.5  at  reduced  velocities  below  and  above  the  flutter  point. 
At  Vr  =  6.2,  below  the  flutter  point,  two  responses  are  possible.  With  an  initial  pitch 
of  a(0)  =  5°,  Figure  2.5(a),  the  airfoil  stabilizes  at  the  stationary  solution.  However, 
if  the  initial  pitch  is  increased  to  ct(0)  =  10°,  Figure  2.5(b),  the  airfoil  exhibits  a  large 
amplitude  LCO.  At  reduced  velocities  above  the  flutter  point,  Figure  2.5(c),  even  a  small 
initial  pitch,  a(0)  =  1°,  can  lead  to  a  large  amplitude  LCO.  A  bifurcation  diagram  with 
this  subcritical  response  is  shown  in  Figure  2.6. 

The  turning  point  in  Figure  2.6  connects  the  subcritical  branch  to  the  stable  pe¬ 
riodic  motion  branch  [32]  and  occurs  at  Vr  ~  5.9.  The  large  amplitude  LCO  branch 
exhibits  a  hysteresis,  i.  e.  the  stationary  solution  bifurcates  to  the  large  amplitude  LCO 
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0  500  1000  1500  2000 


(a)  Vr  =  6.2,  a(0)  =  5° 


(b)  Vr  =  6.2,  a(0)  =  10° 


(c)  Vr  =  6.3,  a(0)  =  1° 


Figure  2.5  Time  histories  of  pitch  with  /3a  =  —3 


Reduced  Velocity 


Figure  2.6  Bifurcation  diagram  of  a  subcritical  response  with  a  turning  point 
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branch  at  the  flutter  point  but  the  large  amplitude  LCO  branch  bifurcates  to  the  sta¬ 
tionary  solution  at  the  turning  point.  This  hysteresis  is  not  simply  a  pure  mathematical 
construct,  the  same  behavior  has  been  observed  in  wind  tunnel  tests  of  a  large  aspect 
ratio  wing  [4]  as  well  as  in  flight  test  of  fighter  aircraft  [2],  The  subcritical  response  with 
a  turning  point  presents  a  greater  risk  to  an  aircraft  than  a  supercritical  response  since 
a  very  large  amplitude  LCO  can  occur  below  the  classically  defined  flutter  point. 

2.5  Uncertainty  Quantification 

While  the  analysis  in  the  previous  section  was  primarily  deterministic,  the  initial 
pitch  angle  a(0)  was  carefully  chosen,  particularly  in  the  case  of  the  subcritical  response, 
so  that  the  desired  behavior  could  be  demonstrated.  The  initial  pitch  angle  plays  a  key 
role  in  the  large  time  behavior  of  this  system.  Thus,  the  initial  pitch  is  chosen  as  a 
stochastic  input. 

It  was  also  demonstrated  in  the  previous  section  that  the  value  of  the  cubic  spring 
constant  played  a  significant  role  in  the  dynamics  of  the  response.  While  the  location 
of  the  flutter  point  does  not  change  with  the  introduction  of  the  nonlinear  structural 
parameters,  the  location  of  the  turning  point  can  change  as  these  structural  parameters 
are  varied.  For  this  reason,  the  second  uncertainty  parameter  chosen  in  this  study  is  the 
cubic  spring  constant,  f3a. 

Before  analyzing  how  these  uncertainties,  a(0)  and  /3a,  propagate  in  time,  a  distri¬ 
bution  of  the  uncertainties  must  be  selected.  While  any  distribution  can  be  chosen  for 
this  model  problem,  the  most  common  starting  point  for  stochastic  projection  methods 
such  as  the  polynomial  chaos  expansion  [15]  is  the  Gaussian  distribution.  The  Gaussian 
distribution  allows  a  random  variable  to  be  expressed  in  terms  of  its  mean  value  and 
its  standard  deviation.  For  example,  realizations  of  the  initial  pitch  angle  and  the  cubic 
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spring  constant  are  obtained  from 


a(0)  =  a(0)  +  £i<5(0) 


(2.14a) 


Pa  ~  Pa  +  £,2$ a 


(2.14b) 


where  the  barred  (“)  quantities  represent  mean  values,  the  tilde  (~)  quantities  represent 
standard  deviations,  and  and  £2  are  zero  mean,  unit  variance  Gaussian  random  vari¬ 
ables.  The  mean  value  in  the  initial  pitch  angle  is  chosen  as  a(0)  =  0.  In  order  to  obtain 
the  desired  responses  shown  in  the  previous  section,  a  large  initial  pitch  angle  is  required. 
This  large  initial  pitch  angle  is  generated  by  selecting  <5(0)  =  0.2  rads  (~  11.5°).  The 
mean  value  of  the  cubic  spring  constant  determines  whether  a  supercritical  or  subcritical 
response  is  produced.  That  is, 


Pa 


{3  generates  a  supercritical  response 
—3  generates  a  subcritical  response 


(2.15) 


The  standard  deviation  of  the  cubic  spring  constant  was  chosen  as  (3a  =  0.3  for  both 
the  supercritical  and  subcritical  responses.  It  should  be  emphasized  that  the  parameters 
were  chosen  only  to  observe  the  desired  responses. 

Monte  Carlo  simulations  of  the  nonlinear  aeroelastic  system  are  now  performed 
in  order  to  quantify  the  effects  of  the  uncertainties.  An  MGS  begins  by  choosing  a 
sample  size,  N,  and  then  obtaining  N  realizations  of  ^  and  £2  based  on  the  Gaussian 
distribution.  These  realizations  of  the  random  variables  provide  realizations  of  a(0)  and 
Pa  from  Eqs.  (2.14a)  and  (2.14b).  Time  integrations  of  the  governing  equations  are  then 
performed  N  times  until  a  fully  developed  response  (stationary  or  LCO)  is  obtained.  The 
responses  are  saved  as  realizations  and,  through  a  Parzen  windowing  approach  [36] ,  used 
to  estimate  the  PDF  of  the  LCO  response. 


2-13 


9  10  11  0  10  20 

aLC0  (degrees)  aLCO  (degrees) 

(a)  pa  =  3,  Vr  =  6.5  (b)  /3a  =  -3,  Vr  =  6.2 

Figure  2.7  Sample  PDFs 


The  Parzen  windowing  approach  provides  a  smooth  estimate  of  the  PDF  through 
the  following  procedure.  The  maximum  and  minimum  values  of  the  LCO  response  are 
first  determined.  A  buffer  value,  A,  is  computed  as  ten  percent  of  the  maximum  value. 
Equally  spaced  nodes  are  selected  on  the  interval  [{ollco) min  —  0.5 A,  (ctlco) max  +  0.5 A], 
where  ollco  is  the  pitch  LCO  response.  Gaussian  distributions  of  a  specified  standard 
deviation,  er,  are  placed  over  each  of  the  nodes  and  the  Gaussian  densities  of  the  responses 
are  averaged  together  [36].  The  standard  deviation  chosen  was  a  =  0.005(az/co)max-  The 
choice  in  the  interval  and  cr  ensured  the  area  under  the  PDF  was  identically  one.  With 
the  smooth  approximation  to  the  PDF,  the  trapezoidal  rule  was  used  to  obtain  the  area 
under  the  PDF.  Sample  PDFs  at  Vr  =  6.5  for  the  supercritical  response  and  Vr  =  6.2  for 
the  subcritical  case  are  shown  in  Figures  2.7(a)  and  (b),  respectively,  ft  should  be  noted 
that  for  the  subcritical  case,  the  finite  width  of  the  PDF  at  cxlco  =  0  is  a  result  of  the 
smoothing  accomplished  by  the  Parzen  windowing  approach. 

Large  time  integrations  were  required  to  ensure  the  responses  obtained  were  fully 
developed.  In  order  to  increase  the  efficiency  of  the  MCS,  i.e.,  decrease  the  required 
computer  run  times,  the  minimum  sample  size  is  selected  that  provides  an  accurate  PDF 
of  the  response.  One  method  for  determining  the  minimum  sample  size  required  is  to 
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run  increasingly  larger  sample  sizes  and  check  for  a  convergence  of  the  mean  value,  /i, 
of  the  response  [37].  An  MCS  was  run  on  the  supercritical  case  with  Vr  =  6.5  and  the 
subcritical  case  with  Vr  =  6.2.  Computer  run  times  and  the  mean  LCO  response  are 
shown  in  Table  2.1.  Each  MCS  for  the  supercritical  case  was  integrated  to  a  maximum 
time  of  rmax  =  3000.  This  was  to  ensure  extremely  small  initial  pitch  angles  achieved 
a  fully  developed  LCO.  For  the  subcritical  response,  a  maximum  time  of  rmax  =  10000 
was  required  to  achieve  either  a  fully  developed  stationary  response  or  a  fully  developed 
LCO.  From  Table  2.1,  it  was  determined  that  N  =  4000  would  be  a  sufficient  sample 
size  to  obtain  accurate  PDFs  of  both  the  supercritical  and  subcritical  responses  with  the 
least  amount  of  computational  effort. 

The  effect  of  the  uncertainties  on  the  bifurcation  behavior  is  demonstrated  for  the 
supercritical  case  in  Figure  2.8(a)  and  (b).  The  supercritical  response  now  has  a  range 
of  possible  LCO  amplitudes  for  each  reduced  velocity  owing  to  the  uncertainty  in  the 
cubic  spring  constant.  The  PDFs  allow  an  estimate  of  failure  for  the  aeroelastic  system. 
Failure  is  defined  as  the  probability  of  encountering  an  LCO.  Due  to  the  smoothing  of  the 
Parzen  windowing  approach  near  cxlco  —  0  for  the  subcritical  case,  failure  is  determined 
from  a  pitch  LCO  amplitude  exceeding  ±1°.  The  area  under  the  PDF  from  oilco  =  1°  to 
the  maximum  LCO  amplitude  yields  the  probability  of  failure.  The  probability  of  failure 
for  the  supercritical  case  is  100%  at  and  above  the  flutter  point. 


Table  2.1  Convergence  of  the  mean  on  MCSs  of  supercritical  and  subcritical  responses 


N 

Supercritical  (rmax  =  3000) 
Run  time  (secs) 

h 

Subcritical  (rmax  =  2000) 
Run  time  (secs) 

R 

1000 

27.98 

9.03 

154.1 

13.9 

2000 

55.86 

9.04 

307.9 

13.2 

3000 

83.88 

9.04 

460.0 

13.4 

4000 

112.1 

9.05 

614.4 

13.5 

5000 

140.1 

9.04 

767.7 

13.6 

6000 

169.8 

9.04 

922.5 

13.6 
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(a)  PDFs  at  various  reduced  ve-  (b)  Bifurcation  diagram 

locities 


Figure  2.8  Effects  of  uncertainties  on  the  supercritical  response 


(a)  PDFs  at  various  reduced  ve¬ 
locities 


(b)  Bifurcation  diagram 


Figure  2.9  Effects  of  uncertainties  on  the  subcritical  response 


The  effects  of  the  uncertainties  on  the  subcritical  response  are  shown  in  Figs.  2.9(a) 
and  (b).  Again,  the  uncertainty  in  the  cubic  spring  constant  leads  to  a  range  of  LCO 
amplitudes.  Introducing  uncertainty  into  the  nonlinear  system  affects  the  location  of 
the  turning  point,  which  now  has  a  probability  of  occurring  as  low  as  Vr  ~  5.7.  This 
variation  in  the  location  of  the  turning  point  agrees  well  with  the  results  reported  by 
Beran  and  Pettit  [8].  With  the  same  definition  of  failure  given  above,  the  probability  of 
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Figure  2.10  Probability  of  failure  for  the  subcritical  response 

failure  for  the  subcritical  case  is  shown  in  Figure  2.10.  While  LCO  below  the  flutter 
point  was  introduced  by  the  nonlinear  restoring  forces,  the  initial  condition  uncertainty 
and  parametric  uncertainty  increase  the  range  where  failure  can  occur  below  the  flutter 
point. 

2. 6  Summary 

In  this  chapter,  a  nonlinear  aeroelastic  system  was  described  that  leads  to  supercrit¬ 
ical  or  subcritical  responses  based  upon  the  structural  parameters  chosen.  Uncertainty 
was  introduced  in  the  analysis  by  assigning  a  Gaussian  distribution  to  the  initial  pitch 
angle  and  the  cubic  spring  constant  in  pitch.  The  response  examined  was  the  LCO  am¬ 
plitude  in  pitch.  It  was  shown  that  the  flutter  point  is  a  function  of  the  linear  system 
only,  while  the  turning  point  of  the  subcritical  response  is  affected  by  both  nonlinearities 
and  uncertainties.  The  turning  point  is  associated  with  LCO,  and  structural  failure, 
at  reduced  velocities  below  the  flutter  point.  Thus,  to  predict  probability  of  failure, 
both  nonlinearities  and  uncertainties  must  be  included  in  the  design  analysis.  The  next 
chapter  presents  a  new  stochastic  algorithm  that  can  efficiently  quantify  the  effects  of 
uncertainties  in  the  nonlinear  aeroelastic  system. 
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III.  Development  and  Implementation  of  the  Stochastic  Projection 

Method  via  B-Splines 

3. 1  Introduction 

The  previous  chapter  presented  an  MCS  analysis  to  determine  the  probability  of 
failure  of  an  aeroelastic  system.  While  the  MCS  approach  is  robust  enough  to  capture 
both  supercritical  and  subcritical  bifurcations,  it  is  computationally  inefficient.  In  this 
chapter,  the  stochastic  projection  of  the  pitch  LCO  response,  ctLco,  is  examined  with 
respect  to  the  stochastic  axes.  This  examination  reveals  that  the  underlying  behavior 
of  a  bifurcation  in  the  stochastic  domain  is  discontinuous  and  not  amenable  to  spectral 
approaches,  such  as  the  PCE  and  FCE  methods.  A  B-spline  algorithm  for  the  stochastic 
projection  method  is  developed  and  applied  to  the  problem  of  the  structurally  nonlin¬ 
ear  airfoil  with  linear  aerodynamics  studied  in  the  previous  chapter.  It  is  shown  that 
this  B-spline  stochastic  algorithm  is  both  efficient  and  robust  enough  to  deal  with  the 
discontinuous  behavior  associated  with  bifurcations. 

3.2  Stochastic  Projection 

A  stochastic  projection  is  a  projection  of  the  response  onto  the  random  variables. 
Typically  this  projection  is  obtained  at  a  specific  instance  in  time.  For  a  time  domain 
analysis,  at  a  very  large  time,  these  projections  can  become  so  complicated  (nonlinear, 
discontinuous)  in  form  that  even  a  very  large  order  expansion  cannot  capture  the  shape  of 
the  projection  [26].  Beran  and  Pettit  [8]  circumvent  this  issue  with  the  non- intrusive  PCE 
approach.  The  non-intrusive  PCE  approach  essentially  removes  the  time  dependency  by 
obtaining  samples  of  the  LCO  response  once  the  LCO  is  fully  developed.  Beran  and 
Pettit’s  [8]  results  were  obtained  from  an  efficient  cyclic  approach,  not  involving  a  time 
domain  analysis.  This  present  work  examines  the  stochastic  projections  from  a  time 
domain  approach. 


3-1 


7  ^  Stochastic  projection 


Singular  point 


(a)  Vr  =  6.5,  pa  =  3.0,  pa  =  0.0, 
a(0)  =  0,  <5(0)  =  0.2 


(b)  Vr  =  6.5,  Pa  =  3.0,  Pa  =  0.3, 
a(0)  =  0.261779,  a(0)  =  0 


Figure  3.1  Stochastic  projection  of  a^co  onto  the  random  variables  £1  and  £2  (super¬ 
critical  response) 


To  better  understand  the  behavior  of  the  stochastic  projections,  two  projections 
are  plotted  in  Figure  3.1.  These  projections  were  obtained  from  a  four  thousand  MCS 
for  the  supercritical  case  {J3a  =  3.0)  at  Vr  =  6.5.  The  responses  were  sorted  with  respect 
to  each  random  variable.  In  Figure  3.1(a)  a  four  thousand  run  MCS  was  performed  with 
the  cubic  spring  constant  held  at  Pa  =  3  while  the  initial  pitch  angle  was  allowed  to  vary. 
The  stochastic  projection  of  the  response,  a^co ,  with  respect  to  the  uncertainty  in  the 
initial  pitch  angle,  appears  to  be  a  smooth  continuous  function,  but  a  singular  point 
exists  at  £1  =  0.  This  singularity  corresponds  to  the  dynamically  unstable  equilibrium 
branch  on  the  bifurcation  diagram  (Figure  2.3).  The  MCS  was  repeated  holding  the 
initial  pitch  angle  fixed  at  a(0)  =  15°,  to  ensure  LCO  developed,  and  allowing  the  cubic 
spring  constant  Pa  to  vary.  Figure  3.1(b)  shows  the  stochastic  projection  of  the  response 
with  respect  to  the  uncertainty  in  the  cubic  spring  constant,  £2  is  a  smooth  continuous 
function. 

The  stochastic  projections  for  the  subcritical  response  pa  =  —  3  at  Vr  =  6.2  are 
shown  in  Figure  3.2.  In  Figure  3.2(a),  the  stochastic  projection  of  the  response  with 
respect  to  the  uncertainty  in  initial  pitch  indicates  two  jumps  between  a  large  LCO 
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(a)  Vr  =  6.2,  f3a  =  -3.0, 
/ 3a  =  0.0,  a(0)  =  0,  a(0)  =  0.2 


(b)  Vr  =  6.2,  (3a  =  -3.0, 

/3a  =  0.3,  a(0)  =  0.261779, 

a(0)  =  0 


Figure  3.2  Stochastic  projection  of  a^co  onf°  the  random  variables  £i  and  £2  (subcrit- 
ical  response) 


amplitude  and  a  stationary  response.  This  discontinuous  behavior  is  a  manifestation  of 
the  subcritical  branch  with  a  turning  point  as  illustrated  in  Figure  2.6.  Above  a  certain 
absolute  value  of  the  initial  pitch  angle,  an  LCO  develops.  Below  this  absolute  value  of 
the  initial  pitch  angle,  a  dynamically  stable  stationary  solution  is  achieved.  In  Figure 
3.2(b),  with  an  initial  pitch  angle  of  a(0)  =  15°  and  f3a  allowed  to  vary,  the  stochastic 
projection  with  respect  £2  is  a  smooth  continuous  function. 


3.3  Stochastic  Projection  Method  via  B- Splines 

The  stochastic  projections  shown  in  the  previous  section  demonstrate  that  while  a 
polynomial  expansion  may  be  appropriate  for  the  supercritical  case  (Figure  3.1),  provided 
the  singular  point  at  £1  =  0  is  not  included  in  the  development  of  the  expansion,  a 
polynomial  expansion  is  inappropriate  for  the  subcritical  case  (Figure  3.2)  due  to  the 
discontinuous  behavior  in  the  projection  with  respect  to  £1.  Indeed,  any  spectral  approach 
would  have  difficulty  in  resolving  the  discontinuous  behavior. 
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The  stochastic  projections  can  be  treated,  however,  as  piecewise  continuous  over  the 
interval  of  interest.  Splines  would  then  be  an  appropriate  choice  for  approximating  the 
stochastic  projections.  An  important  property  of  the  multivariate  B-splines  chosen  for 
this  work  is  that  they  are  a  compact  support  basis;  that  is,  the  influence  of  any  particular 
B-spline  coefficient  extends  over  a  few  intervals  [30].  The  practical  importance  of  this 
property  is  that  oscillations  in  the  vicinity  of  discontinuous  behavior  can  be  avoided 
with  the  proper  choice  in  the  order  of  the  B-spline.  The  order  used  in  this  work  is 
=  k^2  =  2,  which  is  equivalent  to  a  piecewise  linear  interpolation  [30].  The  expansion 
of  the  univariate  B-splines  for  and  £2  are  given  by 

—  ^  ;a,£i  (^1)  (3-1) 

i=  1 

—  ^  ^  ®j&Bj,k£o,xeo  (£2)  (3-2) 

3= 1 

where  x ^  is  the  vector  of  knots  on  the  £1  axis  and  x^2  is  the  vector  of  N^2  knots 
on  the  £2  axis.  The  multivariate  B-spline  is  the  tensor  product  of  the  two  univariate 
B-splines  and  can  be  explicitly  written  as  [30] 

*ei  Nt2 

<*(6,6)  =  EE  (6)5^,^,  (6).  (3.3) 

i= 1  j= 1 

The  coefficient  matrix,  made  up  of  the  elements  ,  is  solved  by  repeated  evaluations 
of  the  univariate  spline  interpolation  problem,  the  details  of  which  are  given  in  de  Boor 

(30]. 

Now  an  efficient  stochastic  algorithm  that  is  robust  enough  to  identify  bifurcations 
can  be  described.  First,  an  interval  is  chosen  in  the  stochastic  domain  over  which  samples 
of  the  LCO  amplitude  will  be  obtained.  Next,  the  location  of  nodes  to  be  sampled  in 
that  interval  is  determined.  The  multivariate  B-spline  that  approximates  the  response 
surface  in  the  stochastic  domain  is  then  determined.  Finally,  an  MCS  is  performed  on 


«(6) 

«(6) 
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this  response  surface  to  estimate  the  PDF  of  the  response.  Each  of  the  steps  of  this 
algorithm  are  now  described  in  detail. 

The  interval  selected  is  based  on  the  distribution  of  random  variables.  For  this 
work,  the  random  variables  have  a  Gaussian  distribution.  Strictly  speaking,  the  interval 
defined  for  a  Gaussian  random  variable  is  the  whole  real  line.  There  is,  however,  a 
practical  limit  to  this  interval.  The  random  variables  control  parametric  uncertainties 
and  these  parameters  must  be  kept  within  specified  bounds  in  order  for  the  equations  of 
motion  to  be  numerically  well  behaved.  In  a  real  world  application  the  mean  and  standard 
deviation  may  be  dictated  by  manufacturing,  instrument  error,  or  physical  phenomena 
such  as  wind  gusts.  For  this  model  problem  the  mean  and  standard  deviation  were 
selected  to  illustrate  supercritical  and  subcritical  responses.  It  was  found  that  for  the 
large  standard  deviation  chosen  for  the  initial  pitch  angle,  d(0)  =  0.2  radians,  the  model 
problem  became  unstable  for  large  values  of  £1  (£1  >  4).  Since 


/  1  \2  ,4  ,4 

W2 nj  J-4  ./< 


e-3 (*?+**)  (!£,  d£2  «  0.99994 


(3.4) 


99.99%  of  all  responses  occur  on  the  interval  [-4,4]  in  both  stochastic  axes.  Thus,  this 
was  the  interval  selected  for  both  the  and  £2  axes. 

With  the  interval  selected,  nodes  along  the  stochastic  axis  must  be  chosen  that  will 
efficiently  capture  the  desired  projection.  In  order  for  the  B-splinc  to  be  valid  over  the 
entire  interval,  nodes  must  be  selected  that  span  the  entire  interval.  However  a  Gaussian 
distribution  of  the  nodes  would  cluster  the  nodes  near  the  mean  value,  leaving  a  large 
gap  between  the  Gaussian  spaced  nodes  and  the  endpoints  of  the  interval.  Nodes  at  ±4 
are  required  to  span  the  interval.  Nodes  at  ±2.5  are  added  to  bridge  the  gap  between  the 
Gaussian  spaced  nodes  and  the  endpoints  of  the  interval.  The  Gaussian  spaced  nodes  are 
determined  by  integration  of  the  Gaussian  PDF,  w(£i)  or  w(£2)  (Eq.  (1.6)  with  d  =  1). 
That  is,  the  nodes  are  determined  in  one  dimension  and  the  same  nodes  are  used  in  the 
other  dimension.  The  determination  of  the  nodes  proceeds  as  follows.  Some  number  of 
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Table  3.1  Gaussian  distributed  nodes 


1  =  2 

I  =  4 

1  =  8 

- 

- 

±0.15731 

- 

±0.31864 

±0.31684 

- 

- 

±0.48878 

±0.67449 

±0.67449 

±0.67449 

- 

- 

±0.88715 

- 

±1.15035 

±1.15035 

- 

- 

±1.53412 

±2.5 

±2.5 

±2.5 

±4.0 

±4.0 

±4.0 

nodes  are  selected  on  the  interval  [0,4],  Let  this  number  of  nodes  be  represented  by  I 
along  the  G  axis  and  J  along  G  axis.  Noting  that 


d£i 


1 

2 


(3.5) 


the  remaining  nodes  are  selected  based  on  equal  intervals  of  the  probabilities.  For  exam¬ 
ple,  with  1  =  2  and  one  node  fixed  at  £1  =  4,  the  second  node  is  determined  from 


d£i 


1 

4 


(3.6) 


where  a  is  the  location  of  the  node.  A  simple  search  leads  to  a  value  of  a  =  0.67449. 
Proceeding  in  a  like  manner  for  /  =  4  and  1  =  8  leads  to  the  node  values  given  in  Table 
3.1.  Due  to  the  symmetry  of  the  nodes  and  the  extra  nodes  at  £1  =  ±2.5,  the  total 
number  of  nodes  along  the  £1  stochastic  axis  is  given  by  2/  ±  2. 

Samples  of  the  response  are  obtained  at  the  selected  values  of  G  and  G-  From  these 
samples,  the  multivariate  B-spline  problem  can  be  solved  for  the  coefficient  matrix, 
in  Eq.  (3.3).  After  these  coefficients  are  determined,  an  MCS  is  performed  on  Eq.  (3.3). 
This  MCS  is  extremely  efficient,  as  is  demonstrated  in  the  next  section. 

Before  preceeding  to  the  next  section,  it  is  noted  that  this  algorithm  fits  into  the 
stochastic  projection  framework  through  the  following  derivation.  Consider  the  set  of 
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nodes  x  such  that  xt  e  [—a,  a]  and  x\  =  —a  <  x2  <■■■<  xj  =  a  for  some  integer 
/.  Since  it  was  noted  that  the  above  algorithm  is  a  piecewise  linear  interpolation,  the 
appropriate  basis  for  expansion  would  be  the  hat  function,  which  has  the  property  [30] 


^  ji.Xi)  $ij 


(3.7) 


where  8tl  is  the  Dirac  delta  function.  Thus 


1  i  —  j 
0  i^j 


(3.8) 


Also  notice  that  8j(x)  =  8(x  —  Xj),  which  indicates  sampling  at  a  node.  Then  inner 
products  with  w(x)  =  1  are  given  as 


(*,,6,)=  [  %(x)6(x~xi)dx  =  %(xi)  =  6ij 


(3.9) 


r 

(xTj,  Sj)  =  /  x^i(x)S(x  —  Xj)  dx  =  Xjd>i(xj)  =  XjSij 

J  —a 

The  piecewise  linear  approximation  to  the  response  a(t,x )  becomes 


(3.10) 


i 

a(t,  x)  =  y,  oiiit^^x)  (3.11) 

1=1 


and  can  be  substituted  into  Eq.  (1.3).  By  performing  the  stochastic  projection  as  de¬ 
scribed  by  Le  Maitre  [18]  and  making  use  of  the  above  inner  products,  the  following 
equation  results 

C[aj(t)]  +  (j3a  +  Xjfio^J  aj(t)  =  0  (3.12) 

Thus,  the  samples  of  the  random  variable  are  used  to  obtain  samples  of  the  response  and 
a  piecewise  linear  approximation  to  the  response  is  obtained  from  the  B-splines.  While 
the  PCE  method  provides  a  best  fit  to  the  stochastic  projection  and  converges,  in  the 
mean  square  sense  to  an  MCS,  [38]  the  piecewise  linear  interpolation  is  exact  at  the 
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nodes  and  converges  exactly  to  the  MCS  when  the  number  of  nodes  equals  the  number 
of  samples  obtained  by  the  MCS. 

3.4  Results 

The  stochastic  algorithm  is  now  applied  to  the  aeroelastic  system  described  in  the 
previous  chapter.  The  supercritical  case  is  examined  first.  Setting  1  =  2  and  holding  the 
cubic  spring  constant  at  (3a  =  3.0,  samples  were  obtained  along  the  £2  axis  at  Vr  =  6.5  to 
determine  the  coefficients  in  Eq.  (3.1).  An  MCS  on  Eq.  (3.1)  produced  the  stochastic 
projection  shown  in  Figure  3.3(a).  Setting  J  =  2  and  holding  the  initial  pitch  constant 
at  a(0)  =  15°,  samples  were  obtained  along  the  £2  axis.  The  univariate  spline,  Eq.  (3.2) 
was  solved  for  the  coefficients  ctj^2  and  an  MCS  on  this  equation  produced  the  stochastic 
projection  shown  in  Figure  3.3(b).  These  projections  are  in  excellent  agreement  with  the 
projections  shown  in  Figure  3.1.  This  should  be  expected  since  the  projections  are  linear 
or  very  nearly  linear  across  the  respective  intervals. 


9 

8 


(a)  Vr  =  6.5,  (3a  =  3.0,  /3a  =  0.0, 

a(0)  =  0,  <S(0)  =  0.2,  I  =  2 


5, 


(b)  Vr  =  6.5,  f3a  =  3.0, 

j3a  =  0.3,  a(0)  =  0.261779, 

a(0)  =  0.0, J  =  2 


Figure  3.3  Stochastic  projection  of  «lco  onto  the  random  variables  £1  and  £2  (super 
critical  response) 
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Figure  3.4  PDF  comparisons  for  the  supercritical  response 

Now  allowing  both  the  initial  pitch  angle  and  the  cubic  spring  constant  to  vary  and 
setting  1  =  2  and  J  =  2,  thirty-six  samples  of  the  LCO  responses  were  obtained.  After 
the  coefficients  in  Eq.  (3.3)  were  determined,  an  MCS  was  performed  on  this  equation. 
The  PDF  of  the  response  was  computed  and  compared  to  the  PDF  from  the  MCS  on 
the  governing  equations.  The  agreement  between  the  two  PDFs  is  excellent,  as  shown  in 
Figure  3.4.  Again  this  should  be  expected  due  to  the  smoothness  of  the  approximated 
functions.  Obtaining  samples  of  the  response,  solving  the  multivariate  B-spline  problem, 
and  performing  an  MCS  on  Eq.  (3.3)  was  two  orders  of  magnitude  faster  than  performing 
an  MCS  on  the  governing  equations. 

The  subcritical  response  was  examined  at  Vr  =  6.2.  The  stochastic  projection  in 
Figure  3.5(a)  was  obtained  as  before  with  the  initial  pitch  angle  allowed  to  vary  and  the 
cubic  spring  constant  held  at  /3a  =  —3.  The  stochastic  projection  in  Figure  3.5(b)  was 
obtained  by  allowing  the  cubic  spring  constant  to  vary  and  keeping  the  initial  pitch  at 
a(0)  =  15°.  The  projection  with  respect  to  required  1  =  8  (18  samples)  to  obtain 
a  good  approximation  to  the  discontinuous  behavior.  The  projection  with  respect  to  £2 
required  J  =  4  (10  samples)  to  obtain  the  proper  curvature.  These  stochastic  projections 
compare  well  with  those  in  Figure  3.2. 

Again  allowing  both  the  initial  pitch  angle  and  the  cubic  spring  constant  to  vary, 
PDFs  were  obtained  from  the  stochastic  projection  method  with  I  =  2,4,  and  8  and 
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(a)  Vr  =  6.2,  f3a  =  -3.0, 
/ 3a  =  0.0,  d(0)  =  0,  a(0)  =  0.2, 
I  =  8 


(b)  Vr  =  6.2,  /3a  =  -3.0, 

/3a  =  0.3,  a(0)  =  0.261779, 

d(0)  =  0.0,  J  =  4 


Figure  3.5  Stochastic  projection  of  onto  the  random  variables  £i  and  ^2  (subcrit- 
ical  response) 


Figure  3.6  PDF  comparisons  for  the  subcritical  response 

J  =  4,  (60,  100,  and  180  samples,  respectively).  These  PDFs  are  compared  with  the  PDF 
from  the  MCS  and  are  shown  in  Figure  3.6.  As  the  number  of  nodes  increase  along  the  /p 
axis,  agreement  with  the  PDF  obtained  from  the  MCS  improves  rapidly.  The  PDF  with 
/  =  4,  is  already  in  excellent  agreement  with  the  MCS  results.  With  I  =  8,  the  steeper 
slope  in  the  approximation  of  the  discontinuous  behavior  results  in  fewer  realizations 
between  the  peak  and  secondary  responses  (a  better  bi-modal  approximation). 
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(a)  Supercritical  response  (b)  Subcritical  response 

Figure  3.7  Bifurcation  diagrams  built  with  the  stochastic  projection  method 

With  the  stochastic  algorithm  verified  for  the  above  cases,  bifurcation  diagrams 
with  the  inclusion  of  uncertainties  were  developed.  For  f3a  =  3,  PDFs  were  obtained  over 
the  same  range  of  reduced  velocities  as  for  the  MCSs.  Setting  I  =  J  =  2  was  sufficient  for 
these  cases  since  all  the  responses  were  supercritical.  For  j3a  =  —3,  the  PDFs  below  the 
flutter  point  were  obtained  with  1  =  8  and  J  =  4,  since  these  responses  were  subcritical. 
Above  the  flutter  point,  once  again  I  =  J  =  2  was  sufficient.  These  PDFs  were  then 
used  to  build  the  bifurcation  diagrams  shown  in  Figure  3.7.  For  the  supercritical  case, 
the  mean  pitch  LCO  amplitudes  from  the  stochastic  projection  method  are  within  2% 
of  the  amplitudes  from  the  MCS  results  (see  Figure  2.8(b)).  The  mean  pitch  LCO 
amplitudes  for  the  subcritical  case  were  much  more  difficult  to  determine.  Due  to  the 
lack  of  resolution  in  the  £i  axis,  numerous  responses  are  realized  between  the  stationary 
response  and  the  maximum  pitch  LCO  amplitude  (see  Figure  3.6).  A  minimum  pitch 
LCO  amplitude  could  not  be  determined  below  Vr  =  5.902.  Thus,  while  the  stochastic 
algorithm  accurately  predicts  an  LCO  response  below  Vr  =  5.902,  only  the  maximum 
value  of  that  response  could  be  determined.  Above  Vr  =  5.902,  the  mean  pitch  LCO 
amplitudes  were  within  2%  of  the  amplitudes  predicted  by  the  MCS  analysis  (see  Figure 
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Figure  3.8  Comparison  of  probability  of  failure  predictions  with  MCS  and  B-spline 
algorithm 

2.9(b)).  It  should  be  noted  that  the  time  required  to  build  these  bifurcation  diagrams 
was  approximately  two  orders  of  magnitude  less  than  with  the  MCS  approach. 

From  the  subcritical  responses,  an  estimate  of  the  probability  of  failure  was  ob¬ 
tained.  As  in  the  MCS  analysis,  failure  is  defined  as  encountering  an  LCO.  The  proba¬ 
bility  of  failure  is  computed  by  integrating  the  area  under  the  PDF  from  «lco  =  1°  to 
the  maximum  LCO  pitch  amplitude.  A  comparison  between  the  MCS  prediction  and  the 
B-spline  prediction  is  shown  in  Figure  3.8.  The  presence  of  realizations  between  the  peak 
and  secondary  responses  lead  to  a  conservative  estimate  in  the  probability  of  failure,  as 
seen  in  Table  3.2.  Each  failure  estimate  from  the  MCS  required  over  ten  minutes  of  com¬ 
puter  time,  while  each  estimate  from  the  stochastic  algorithm  required  approximately 
fifteen  seconds. 

3. 5  Summary 

An  efficient  algorithm  for  determining  the  propagation  of  uncertainties  in  a  highly 
nonlinear  aeroelastic  system  has  been  presented.  The  uncertainty  in  the  initial  pitch 
angle  a(0)  was  allowed  to  propagate  in  a  time  dependent  manner  until  an  LCO  or  a 
dynamically  stable  solution  was  achieved.  The  stochastic  algorithm  consisted  of  building 
an  interpolating  function  in  the  stochastic  domain  through  sampling  and  the  use  of  a 
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Table  3.2  Estimates  for  the  probability  of  failure 


Reduced  Velocity  (Vr) 

MCS  (%) 

Stochastic  Algorithm  (%) 

Difference  (%) 

5.70 

0.005 

0.91 

0.905 

5.75 

0.33 

1.02 

0.69 

5.80 

0.70 

2.52 

1.82 

5.85 

3.23 

5.17 

1.94 

5.90 

7.23 

10.9 

3.67 

5.95 

12.4 

17.1 

4.7 

6.00 

19.2 

27.6 

8.4 

6.05 

26.1 

33.8 

7.7 

6.10 

35.5 

38.2 

2.7 

6.15 

45.0 

50.7 

5.7 

6.20 

57.1 

61.4 

4.3 

6.25 

74.1 

77.5 

3.4 

6.279 

100.0 

100.0 

0 

multivariate  B-spline.  Advantages  of  this  method  over  the  traditional  PCE  approach 
were  that  expected  values  were  never  computed  or  stored  and  two  orders  of  magnitude 
reduction  in  computational  time  over  an  MCS  analysis  were  obtained.  Based  on  the  PDFs 
predicted  by  the  stochastic  algorithm,  a  rapid  and  accurate  estimate  of  the  probability 
of  failure  was  also  obtained.  It  should  be  noted  that  while,  thus  far,  only  the  pitch  LCO 
amplitude  has  been  examined  as  a  response,  the  stochastic  algorithm  is  easily  extended 
to  examining  other  responses  such  as  plunge  LCO  amplitude  or  the  coupled  pitch  and 
plunge  frequency.  These  additional  responses  will  be  examined  in  the  next  chapter  where 
a  computation  of  the  nonlinear  aerodynamics  is  obtained  from  the  Euler  equations. 
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IV.  The  Inviscid  Aeroelastic  Code  EULER-AE 

4-1  Introduction 

To  this  point  aerodynamic  nonlinearities  in  the  aeroelastic  system  have  been  ne¬ 
glected.  Limit  cycle  oscillations  often  occur  in  transonic  flow  [2] ,  where  the  nonlinearities 
dne  to  moving  shocks  cannot  be  ignored.  Other  effects,  such  as  boundary  layer  separation 
and  turbulence,  are  also  important  in  the  study  of  LCO  at  transonic  speeds  [5] ;  however, 
many  researchers  have  demonstrated  LCOs  on  typical  airfoil  sections  with  the  inviscid 
Euler  equations  [34,  39,  40,  41].  Therefore,  validation  of  the  stochastic  algorithm  pre¬ 
sented  in  the  previous  chapter  does  not  require  viscosity  to  be  included  in  the  modeling 
of  the  fluid  dynamics. 

An  inviscid  aeroelastic  code,  EULER-AE,  has  been  developed  to  investigate  the 
response  of  a  nonlinear  airfoil  in  an  inviscid  transonic  flow.  The  code  implements  a  finite 
volume  formulation  and  loose  coupling  between  the  structure  and  fluid,  i.e.,  the  structural 
dynamics  lag  the  fluid  dynamics  by  the  integration  time  step  size.  The  airfoil  selected 
for  this  study  was  the  symmetric  NACA  64A006  airfoil.  This  airfoil  was  investigated  by 
Morton  and  Beran  [34]  to  determine  the  flutter  point  of  a  linear  airfoil  (/3a  =  ya  =  0)  in 
inviscid  flow.  The  grid  generation  and  dynamic  boundary  conditions  used  by  EULER-AE 
are  presented.  The  nonlinear  airfoil  equations  are  modified  to  use  the  notation  of  Morton 
and  Beran  [34] .  The  strategy  of  EULER-AE  is  to  first  attempt  convergence  to  a  steady 
state  solution  at  a  fixed  initial  pitch  angle.  Once  convergence  or  the  maximum  number 
of  iterations  is  achieved,  the  airfoil  is  allowed  to  move  subject  to  the  lift  and  moment 
induced  by  the  fluid  flow  (See  Figure  4.1).  Static  solutions  with  the  airfoil  at  different 
angles  of  attack  were  computed  with  EULER-AE  and  compared  against  results  from 
the  commercial  fluid  solver  FLUENT™  as  well  as  experimental  data  [42] .  A  comparison 
with  the  parameters  used  by  Morton  and  Beran  to  obtain  an  LCO  for  this  airfoil  is  also 
presented. 
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Figure  4.1  EULER-AE  Flow  Chart 


4-2  Euler  Equations 

Before  presenting  the  Euler  equations,  a  new  set  of  nondimensional  parameters  are 
introduced.  In  Sec.  2.2  the  nondimensionalization  was  based  on  the  airfoil  midchord,  b, 
as  a  reference  length.  In  CFD  algorithms,  the  airfoil  chord,  c,  is  typically  chosen  as  a 
reference  length.  The  nondimensional  parameters  are  now 
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(4.1) 


From  this  point  the  asterisk  is  dropped  from  the  nondimensional  parameters  and  the 
primitive  variables  and  grid  locations  are  assumed  to  be  nondimensional  unless  otherwise 
stated. 
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The  Euler  equations  in  generalized  coordinates  with  a  moving  grid  are  given  by  [43] 


d_ 

dt 


d 


(4.2) 


where  the  Jacobian  J  is  equivalent  to  the  cell  volume,  the  conserved  variables  are 


Q 


p 

pu 

pv 

E, 


(4.3) 


and  the  flux  vectors  are  given  by 


pU 

pV 

ixP  +  pUu 

,  F  = 

rjxP  +  pVu 

CyP+pUv 

TJyP  +  pV  V 

(p  +  Et)U  -£tp 

(p  +  Et)V  -  ptP 

The  contravarient  velocities  U  and  V  are  obtained  from 


(4.4) 


U  T  t \x U  T  ^yV 


(4.5a) 


V  =  pt  +  T)XU  +  T)  yV 


(4.5b) 


For  a  finite  volume  approach,  the  spatial  metrics  of  transformation  (£x,  £y,  rjx ,  and  rjy) 
are  obtained  from  a  relationship  with  the  area  vectors  normal  to  the  cell  faces  [44],  The 
temporal  metrics  of  transformation  are  obtained  from  a  first  order  backward  difference. 
At  time  level  n,  the  x-component  and  ^/-component  of  the  grid  velocity  are  given  by 


^.n  ^.n— 1 


Xt  = 


Vt  = 


At 

y1^ _ y^ 


At 


(4.6a) 

(4.6b) 
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respectively,  and  the  metrics  are  computed  from 


6  =  -  {£xXt  +  ZyVt)  (4.7a) 

Vt  =  ~  (VxXt  +  VyVt)  (4.7b) 


An  equation  of  state  is  required  to  close  the  system  of  equations.  For  an  ideal  gas 
with  a  constant  ratio  of  specific  heats,  7,  the  equation  of  state  is  written  as 

Et  =  — ( u 2  +  v 2)  (4.8) 

7  —  1  2  v  7 

where  7  =  1.4. 


4-3  Numerical  Procedure 


The  two-stage  Runge-Kutta  time  integration  in  EULER-AE  uses  Euler’s  forward 
difference  as  a  predictor  and  the  trapezoidal  method  as  a  corrector.  Atkinson  [29]  shows 
this  is  scheme  is  second  order  accurate  in  time.  Notationally,  the  integration  scheme  can 
be  written  as 


Q?y  =  Q“i  + 


1  At 


hO 


hJ 


2  V, 


(Ki  +  Ki) 


(4.9) 


where  Vjj  is  the  cell  volume,  7 Zni:j  is  the  residual  at  time  step  n,  and  7?.\j  is  the  residual 
obtained  from  the  Euler  prediction. 


The  residual  IZij  is  computed  from  the  flux  differencing  Roe  method  [45],  a  first 
order  (spatial)  approximate  Riemann  solver.  The  Roe  method  linearizes  the  flux  Jacobian 
matrices  [A]  and  [B],  where 


<9E 

dQ 


(4.10) 


by  replacing  these  matrices  with  Roe  averaged  matrices  [A]  and  [B]  so  that  the  solution 
to  the  linear  problem  provides  an  approximate  solution  to  the  nonlinear  problem  [44], 
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The  averaging  is  based  on  the  states  to  the  left  and  right  of  the  node  being  computed, 
e.g. 


Q  L  =  Qi-1  (4.11a) 

Q  r  =  Q  i  (4.11b) 

along  a  constant  j  line.  The  Roe  averaged  matrices  are  represented  by  the  diagonal 
transformations 


M  =  RHIAellPr1  (4.12a) 

[B]  =  [r,][|A,|][r,]_1  (4.12b) 


where  [T^]  is  composed  of  columns  of  the  right  eigenvectors  of  [A],  [ T, t]_1  is  composed  of 
rows  of  the  left  eigenvectors  of  [A],  and  [|A^|]  is  a  diagonal  matrix  of  the  absolute  values 
of  eigenvalues  of  [A].  The  diagonalization  for  [B\  is  made  up  of  similar  matrices.  An 
efficient  formulation  of  the  Roe  method  that  pre-multiplies  the  transformation  matrices 
[45]  is  used  in  EULER- AE. 

The  residuals  are  computed  by  summing  the  fluxes  across  the  cell  faces  from 


where 


7? 


hi 


+  F, , 


ti+5 


(4.13) 


Fr  +  Fl  -  [^[[A.O^-^Qr  -  Ql) 


(4.14a) 

(4.14b) 


Higher  order  spatial  discretizations  can  be  obtained  from  the  monotone  upstream- 
centered  schemes  for  strong  conservation  laws  (MUSCL)  [44],  The  MLISCL  scheme  ex- 
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trapolates  left  and  right  states  from  more  than  one  neighboring  node.  For  example 


Ql  =  Q,-,  +  LjAQi-i  + 

(4.15a) 

Q r  —  Qi  ^  <5Qi_§  ^  5Qj+i 

(4.15b) 

where  6  is  the  central  difference  operator  [44],  The  constant  k  determines  the  order  and 
method  of  the  MUSCL  scheme  from 


k  = 


-1 

0 


l 

3 


1 


2nd  order  upwind  scheme 
2nd  order  Fromm’s  method 
3rd  orc}er  upwind  scheme 
2nd  order  central  difference 


(4.16) 


The  third  order  Roe  method  was  used  throughout  this  investigation. 


Because  of  the  extrapolations  made  to  obtain  higher  order  schemes  non-physical 
overshoots,  such  as  negative  densities,  are  numerically  possible  in  regions  with  discon¬ 
tinuities,  i.e.  shocks.  Limiters  are  used  to  prevent  these  overshoots  by  enforcing  the 
first  order  Roe  method  in  the  vicinity  of  discontinuities.  ELILER-AE  uses  the  MINMOD 
limiter  [44]  to  control  overshoots. 


For  steady  state  convergence,  local  time  stepping  is  used  as  an  acceleration  tech¬ 
nique.  The  time  step  A tij  is  determined  from  each  cell  volume  and  local  contravarient 
velocities  from 


a.  CFL 

,J  ~  ‘T[(v)«  +  (a2)„] 


(4.17) 


where 


^ i,j 

(-^2 )i,j  ~  ui,j 


(4.18a) 

(4.18b) 
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and  CFL,  the  Courant-Friedrichs-Lewey  number,  is  a  multipication  factor  that  is  chosen 
to  achieve  the  largest  possible  explicit  time  step  without  losing  stability.  For  the  third 
order  Roe  scheme  employed,  CFL  =  1  was  used. 

A  time  accurate  integration  requires  a  global  time  step.  This  global  time  step  is 
the  minimum  time  step  obtained  from  Eq.  (4.17).  For  this  study,  a  constant  time  step 
of  At  =  0.0002  was  chosen.  This  was  equivalent  to  using  CFL  «  1.25  in  Eq.  (4.17). 

4-4  Grid  Generation 

Gordnier  and  Melville  [7]  demonstrated  the  ability  to  predict  LCO  behavior  is 
dependent  upon  the  fidelity  of  the  grid.  They  showed  that  for  an  inviscid  solver  the 
estimated  flutter  point  was  in  better  agreement  with  experimental  data  with  a  coarse 
grid  rather  than  with  a  fine  grid.  Since  boundary  layer  effects,  such  as  flow  separation 
and  turbulence,  cannot  be  modeled  by  an  inviscid  code,  the  shock  is  predicted  in  the 
wrong  location.  A  high  order  spatial  integration  on  a  fine  grid  results  in  the  shock  being 
sharply  refined  at  this  wrong  location,  leading  to  errors  in  the  lift  and  moment  coefficients. 
While  grid  refinement  was  an  issue  near  the  leading  edge  of  the  NACA  64A006  airfoil  due 
to  the  small  leading  edge  radius,  the  grid  was  allowed  to  be  coarser  towards  the  trailing 
edge  in  order  to  spread  (or  weaken)  the  shock  over  a  larger  area. 

The  ordinates  for  the  NACA  64A006  airfoil  are  shown  in  Table  4.1  [46,  47].  In  order 
to  construct  the  airfoil,  a  cubic  spline  of  the  ordinates  was  merged  to  the  leading  edge 
radius.  The  trailing  edge  was  pinched  down  so  that  the  trailing  edge  radius  was  zero. 
This  was  necessary  since  the  inviscid  code  cannot  handle  the  flow  separation  induced  by 
the  small  trailing  edge  radius.  The  change  in  geometry  is  minor  and  should  not  have  a 
significant  effect  on  the  aerodynamics. 

A  257  x  65  point  C-gricl  was  constructed  from  the  airfoil  geometry.  A  cut  line  was 
defined  from  the  airfoil  trailing  edge  to  the  the  far  held  boundary.  Along  the  j  —  1  line, 
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Tabic  4.1  Ordinates  for  the  NACA  64A006  Airfoil3, 


X 

y 

X 

y 

0.0000 

0.00000 

0.2000 

0.02557 

0.0005 

0.00158 

0.2500 

0.02757 

0.0010 

0.00221 

0.3000 

0.02896 

0.0015 

0.00270 

0.3500 

0.02977 

0.0020 

0.00311 

0.4000 

0.02999 

0.0025 

0.00347 

0.4500 

0.02945 

0.0030 

0.00379 

0.5000 

0.02825 

0.0035 

0.00409 

0.5500 

0.02653 

0.0040 

0.00435 

0.6000 

0.02438 

0.0050 

0.00485 

0.6500 

0.02188 

0.0075 

0.00585 

0.7000 

0.01907 

0.0125 

0.00739 

0.7500 

0.01602 

0.0250 

0.01016 

0.8000 

0.01285 

0.0500 

0.01399 

0.8500 

0.00967 

0.0750 

0.01684 

0.9000 

0.00649 

0.1000 

0.01919 

0.9500 

0.00331 

0.1500 

0.02283 

1.0000 

0.00013b 

aL.E.  radius  0.00246,  T.E.  radius  0.00014 
bThis  y-ordinate  was  changed  to  0.00000 


the  cut  line  and  airfoil  surface  are  given  by 


i 


1-65 
65  -  129 

< 

129  -  193 
193  -  257 


cut  line 

airfoil  lower  surface 
airfoil  upper  surface 
cut  line 


(4.19) 


The  leading  edge  radius  required  fine  spacing  to  reduce  pressure  oscillations  seen  in  a 
coarser  grid.  The  initial  spacing  on  the  leading  edge  was  0.00025  chords  and  the  spacing 
was  allowed  to  increase  slowly  as  the  curvature  decreased.  The  grid  spacing  was  coarsest 
near  the  midpoint  of  the  airfoil  and  then  became  more  finely  spaced  near  the  trailing 
edge  to  accommodate  flow  field  gradients  in  this  region.  The  final  grid  point  spacing 
at  the  trailing  edge  was  approximately  0.01  chords.  The  radial  spacing  was  chosen 
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Figure  4.2  NACA64A006  structured  grid  -  far  field 


(a)  Leading  edge  (b)  Trailing  edge 

Figure  4.3  NACA64A006  structured  grid  -  detail 
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to  minimize  cell  aspect  ratio  near  the  leading  edge  without  becoming  so  fine  as  to  be 
considered  a  viscous  mesh.  The  initial  grid  spacing  off  the  wall  was  0.001  chords  and 
spacing  exponentially  increased  out  to  fifteen  chord  lengths. 

Two  control  points  were  defined  approximately  10°  forward  and  aft  of  the  airfoil 
on  the  far  field  boundary.  Grid  points  were  redistributed  between  these  control  points 
until  the  algebraic  grid  was  nearly  orthogonal  to  the  airfoil  surface.  An  elliptic  solver 
was  then  run  on  the  subdomain  defined  by  the  first  forty  grid  points  off  the  wall  and  the 
cut  line  to  further  improve  orthogonality.  The  resulting  grid  is  shown  in  Figure  4.2.  A 
close-up  of  the  grid  near  the  leading  and  trailing  edges  is  shown  in  Figure  4.3. 

4-5  Boundary  Conditions 

The  boundary  conditions  consist  of  three  types  -  solid  wall  (impermeable  surface), 
far  field  (inflow  or  outflow),  and  a  cut.  The  following  first  order  formulation  comes 
primarily  from  Whitfield  and  Janus  [48]  with  modifications  to  include  the  fluxes  induced 
by  the  rigid  body  motion  of  the  grid.  Since  the  flow  solver  is  node  centered,  the  metrics 
of  transformation  at  the  nodes  are  computed  with  standard  finite  difference  techniques 
(see  Tannehill  et  al.  [44]). 

The  inviscid,  solid  wall  boundary  condition  requires  the  flow  to  remain  tangential 
to  the  wall.  For  a  general  airfoil  geometry,  the  tangential  flow  is  obtained  by  computing 
the  unit  normals  at  each  node  from  the  metrics  of  transformation,  i.e. 

(4.20a) 
(4.20b) 


The  tangential  velocities  at  the  solid  wall  are  then  given  by 

^ wall  ^  Uhx 
V wall  —  V  U'fly  . 


nr  = 


n„  = 


r)x 


y/tfZ  +  Vy 
Vy 

y/v'Z  +  Vy 
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(4.21a) 

(4.21b) 


Both  the  pressure  and  density  gradients  were  assumed  to  be  zero  at  the  wall  [43].  Thus 
Pwaii  and  pwaii  were  set  to  the  same  values  as  those  one  grid  point  off  the  wall.  The  total 
energy  at  the  wall  was  computed  from  the  ideal  gas  equation  of  state,  Eq.  (4.8). 

The  subsonic  far  held  boundary  conditions  depends  upon  whether  the  how  is  in  to 
or  out  of  the  computational  domain.  The  standard  unit  normal  direction  is  defined  to 
point  into  the  computational  domain.  This  direction  is  in  the  same  direction  (sign  =  1) 
as  the  unit  normal  determined  from  the  metrics  of  transformation  for  the  i  =  1  and 
j  —  1  faces,  and  in  the  opposite  direction  (sign  =  —  1)  for  the  i  =  imax  and  j  =  jmax 
faces.  Thus,  when  sign  x  U  >  0,  the  how  is  in  to  the  computational  domain,  and  when 
sign  x  U  <  0  the  how  is  out  of  the  computational  domain. 

For  a  subsonic  inflow,  three  of  the  characteristic  variables  are  determined  from  the 
free  stream  conditions,  but  the  fourth  must  be  determined  from  the  computational  do¬ 
main.  For  a  subsonic  outflow,  two  of  the  characteristic  variables  must  be  determined  from 
the  computational  domain  while  the  third  is  determined  from  the  free  stream  conditions 
[44],  Whitfield  and  Janus  [48]  relate  these  characteristic  variables  back  to  the  primitive 
variables.  For  subsonic  inflow  the  boundary  conditions  are  given  by 


Pin 


Pin  Poo 

^  (Poo  +  Pref  +  sign  pref  aref  {u  X  [^OO  ('i'Uef  •£<)]  T  Tly  [  t’oo 

Poo  Pref 

Min  =  Moo  +  sign  nx - 

Pref®  ref 

Poo  Pref 

Min  =  Moo  +  sign  ny - 

Pref  ®  ref 


(4.22a) 

(Mref-Pt)]})  (4.22b) 
(4.22c) 

(4.22c!) 


where  the  reference  values  (pref,  pre f,  etc.)  are  obtained  from  the  hrst  interior  node.  The 
subsonic  outflow  boundary  conditions  are  given  by 


Pout  Poo 


(4.23a) 
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_  Poo  Pref 

Pout  Pref  i  2 

aref 


^out  =  ^ref  +  Sign  nx 


Vout  =  vref  ~h  sign  n, 


Poo  Pref 
c 

Pref  ^ref 
Poo  Pref 

i 

Pref^ref 


(4.23b) 

(4.23c) 

(4.23d) 


In  both  cases,  the  total  energy  was  computed  from  the  ideal  gas  equation  of  state. 


The  boundary  condition  employed  along  the  cut  line  of  the  C-grid  was  a  simple 


average.  For  example 


PiJ  =  2^2j+Pw-ij) 


(4.24) 


Lift  and  Moment  Computations 

The  lift  and  moment  coefficients,  assuming  unit  span,  are  given  by 


Ci  = 


Cm  — 


P  d Ay 


p  (. X  d Ay  -  y  d Ax) 


(4.25a) 

(4.25b) 


These  equations  were  derived  assuming  the  grid  underwent  rigid  body  motion  and  the 
moment  was  taken  about  the  origin.  The  minus  sign  in  the  coefficient  of  lift  equation 
ensures  that  lift  is  positive  up.  The  moment  coefficient  is  measured  positive  nose  up. 
The  components  of  the  area  vector  normal  to  the  airfoil  surface  are  given  by 


cM*  =  dL 
J 

d  4  -  rJy 

dAy  -  j  . 


(4.26a) 

(4.26b) 


By  a  straightforward  application  of  the  trapezoidal  rule,  the  lift  and  moment  coefficients 
were  computed  from 


G  =  -pici.i(^).ia-2£  Pi, 


i=ic  i+l 


-2V  ..  -pw 


j/i,  i 


J  '  *c2jl 


(4.27a) 
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Figure  4.4  Notation  for  the  Pitch  and  Plunge  Airfoil  used  in  EULER-AE 


and 


Cm  Pic  i,l 


V  y\  (  V: 

*ici,l  It)  Uic  1,1  (  t  / 
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*c2— 1 
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Py  \  (Vx 
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(Hx 

Xic  2,1  It  _  ^C2,l  T  , 
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(4.27b) 


The  indices  fci  and  ic2  are  located  on  the  trailing  edge  of  the  airfoil,  where  ici  =  65  and 
ic2  =  193. 


Ah  Aeroelastic  Equations  of  Motion 

The  equations  of  motion  for  the  airfoil  differ  slightly  from  those  presented  in  Sec.  2.2. 
The  nomenclature  shown  in  Figure  4.4  was  used  by  Morton  and  Beran  [34]  and  was  also 
used  in  the  development  of  EULER-AE.  Viscous  dampers  are  now  included  in  both  the 
pitch  and  plunge  axes.  They  are  modeled  by  [34] 

Da  =  2Q0rnrfl)1uj0  (4.28a) 

Dh  =  2( ymuh  (4.28b) 
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The  location  of  the  center  of  gravity  is  measured  from  the  airfoil  leading  edge  and  is 
scaled  by  the  airfoil  chord.  Plunge  is  measured  positive  up  and  is  also  scaled  by  the 
airfoil  chord.  Additionally,  the  nondimensionalization  for  reduced  velocity  is  given  by 


Vr 


14c 


U! aC 


(4.29) 


and  the  time  t  is  given  in  Eq.  (4.1).  Cubic  and  quintic  nonlinear  responses  are  still 
assumed  for  the  pitch  spring.  With  the  above  assumptions  and  the  remaining  nondimen- 
sional  parameters  given  in  Sec.  2.2,  the  equations  of  motion  become 


Xa  ..  .  c 0r  . 

y+—a  +  KyyV  +  4 


Ur 

Vr 


y  = 


/is7T 


-Cl 


x. 


y  +  +  2 (a  +  f3aa3  +  7 Qo;5)  —  —— ~Cm  . 

V  r 


142 


/is7T 


(4.30a) 

(4.30b) 


By  letting 


S  = 


C4 

\S,J 


fa\ 


a 

y 

\y ) 


the  equations  of  motion  can  be  written  in  the  compact  notation 


(4.31) 


[M]S  +  [A'i]S  +  [/i3]S3  +  [/iJS5  =  G 


where  the  mass  matrix  is 


[M\ 


10  0  0 
0^01 
0  0  10 
0  0  xa 


(4.32) 


(4.33a) 
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the  linear  stiffness  matrix  is 


[Ki]  = 


0 

0 

0 
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-1 

0 

0 
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0 
0 


0 

-1 

0 


the  nonlinear  stiffness  matrices  are 
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0 
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0 
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0 
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and  the  aerodynamic  forcing  term  is 


G  = 


/o\ 

/ 

0 
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0 
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4 
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Solving  for  S  yields 


S  =  [MY1  (G  -  [K\]S  -  [Ji3]S3  -  [Ji5]S5) 


(4.33b) 


(4.33c) 


(4.33d) 


(4.34) 
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which  fully  expanded  becomes 


where 


A  A 

82  1 

S2 

rZ-xl  [9l  xafl] 

S3 

S3 

\s4) 

\  r2lx2  [rlfi  xagx]  y 

(4.35) 


/i  =  /  -  4 


Ur 

Vr 


A  /■  <X>r 

S^-KyySi 


9l  ~  9  ~  ^y2  (Sl  Pas\  +  7aSl)  —  7,(a-^j-S2  ■ 


Vr 


(4.36a) 

(4.36b) 


Equation  (4.35)  is  integrated  explicitly  in  time  by  the  two-stage  Runge-Kutta  method 
already  described. 


4-7  Rigid  Body  Motion 

At  each  step  in  time,  the  grid  is  translated  and  rotated  with  respect  to  the  elastic 
axis,  where  the  elastic  axis  is  initially  placed  at  the  origin.  Defining 

A a  =  an-  a11-1  (4.37a) 

A  y  =  yn-yn~1  (4.37b) 


the  translation  is  accomplished  by 


and  the  rotation  by 


V  =9 


n— 1 


+  A  y 


cos(Aa)  sin(Aa) 
-  sin(Aa)  cos(Ao;) 


x 


n—1 


(4.38) 


(4.39) 
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The  transformation  metrics  are  unaffected  by  the  translation,  but  must  be  rotated  by 
the  same  rotation  matrix  in  Eq.  (4.39)  in  order  for  the  area  vectors  to  remain  normal  to 
the  cell  faces. 


4-8  Sample  Cases 

Two  static  cases  and  an  LCO  case  are  presented  as  validation  of  the  research  code. 
The  static  cases  are  compared  with  results  from  FLUENT™  and  with  experimental  data 
[42],  The  first  case  was  run  at  =  0.86  and  a  =  0  and  the  second  case  was  run  at  Mm  = 
0.85  and  a  =  4°.  The  resulting  pressure  distributions  are  shown  in  Figure  4.5.  Both 
EULER-AE  and  FLUENT™  are  in  very  good  agreement  with  one  another.  The  predicted 
shock  locations  are  one  cell  different  between  the  two  codes.  This  is  probably  due  to 
EULER-AE  being  run  with  a  third  order  spacial  scheme,  while  FLUENT™  is  limited  to 
a  second  order  scheme.  The  difference  in  order  of  spatial  discretization  also  explains  why 
the  shock  is  more  sharply  defined  with  EULER-AE.  As  expected,  the  inviscid  shocks 
are  stronger  and  further  aft  on  the  wing  than  the  experimental  data,  especially  at  large 
angles  of  attack  (Figure  4.5(b)).  This  will  lead  to  a  prediction  of  a  much  larger  lift 


(a)  Moo  =  0.86,  a  =  0 


(b)  Moo  =  0.85,  a  =  4° 


Figure  4.5  Pressure  distributions  on  the  NACA  64A006  airfoil 
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coefficient  for  the  inviscid  cases.  While  comparison  with  the  experimental  data  is  only 
fair,  due  to  the  lack  of  viscous  and  turbulence  modeling,  the  code  to  code  comparison 
with  FLUENT™  indicates  that  EULER-AE  properly  predicts  an  inviscid  flow  over  the 
airfoil. 

The  time  accurate  fluid-structure  interaction  was  verified  by  comparing  with  an 
LCO  prediction  given  by  Morton  and  Beran  [34],  A  complete  validation  was  not  possible 
due  to  the  different  formulations  of  ELILER-AE  and  Morton  and  Beran  [34],  EULER- 
AE  employs  a  matched  point  analysis  whereas  the  method  of  Morton  and  Beran  [34] 
employs  a  non-matched  point  analysis.  A  matched  point  analysis  ensures  the  reduced 
velocity,  freestream  air  density  and  freestream  Mach  number  are  consistent  with  standard 
atmospheric  conditions  [49].  A  given  density  and  Mach  number  implies  a  freestream 
velocity  which  in  turn  implies  the  values  of  the  viscous  dampers  and  spring  constants, 
as  shown  in  Eqs.  (4.29)  and  (4.30).  Thus,  the  amplitude  and  frequency  of  the  coupled 
response  from  a  matched  point  and  non-matched  point  analysis  will  not  necessarily  be 
the  same  (see  Appendix  A). 

Morton  and  Beran  [34]  use  the  following  parameters  to  obtain  an  LCO  from  the 
NACA  64A006  airfoil  with  linear  pitch  (fja  =  =  0)  and  plunge  springs. 

xcg  =  0.375  =  -0.25  r2a  =  0.25 

Moo  =  0.85  Us  =  125  ur  =  0.2  (4.40) 

C y  =  0.5  Ca  =  0.5 

At  a  reduced  velocity  of  Vr  —  11  and  standard  sea  level  conditions  {p^  =  101325  N/m2 
and  Poo  =  1.225  kg/m3),  ELILER-AE  predicts  the  pitch  and  plunge  LCOs  shown  in 
Figure  4.6.  Table  4.2  shows  the  differences  between  the  results  from  ELILER-AE  and 
Morton  and  Beran  [34],  The  shorter  period  from  ELILER-AE  indicates  that  the  pitch 
spring  stiffness  for  the  aeroelastic  system  is  greater  than  that  from  Morton  and  Beran 
[34],  However,  the  larger  amplitudes  in  pitch  and,  especially,  plunge  indicate  just  the 
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Figure  4.6  LCO  at  Vr  —  11.0 


Table  4.2  Differences  between  LCOs  from  EULER-AE  and  Morton  and  Beran  [34] 


ollco 

Vlco 

Period 

Flutter  Point 

EULER-AE 

15.8° 

0.176 

32.9 

6.78 

Reference  [34] 

13.0° 

0.020 

61.6 

10.28 

opposite.  The  reason  for  the  apparent  discrepancy,  as  will  demonstrated  in  the  next 
chapter,  is  the  flutter  point  for  the  matched  point  aeroelastic  system  is  located  at  a  much 
lower  reduced  velocity  (Vr  ~  6.78)  than  that  of  the  non- matched  point  aeroelastic  system 
(Vr  =  10.28). 

4-9  Summary 

An  inviscid  aeroelastic  code  has  been  developed  to  investigate  the  response  of  a 
nonlinear  airfoil.  The  equations  and  numerical  procedures  used  by  the  code  EULER- 
AE  were  described.  Sample  solutions  to  verify  and  validate  the  numerical  procedures 
were  presented.  In  the  next  chapter,  this  code  is  used  to  determine  the  effects  of  initial 
condition  and  parametric  uncertainty  on  the  nonlinear  NACA  64A006  airfoil. 
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V.  A  CFD  Application  of  the  Stochastic  Projection  Method  via  Id- 

Splines 

5. 1  Introduction 

A  typical  aeroelastic  simulation  performed  by  EULER-AE  took  approximately  36 
to  48  hours  to  reach  a  fully  developed  LCO,  making  the  MCS  approach  to  determining 
PDFs  of  a  response  impractical.  The  stochastic  algorithm  developed  in  Chapter  111 
was  applied  to  this  CFD  application  to  determine  the  subcritical  response  of  the  NACA 
64A006  airfoil.  The  subcritical  case  was  examined  since  LCOs  can  be  obtained  below 
the  classically  defined  flutter  point. 

To  produce  a  subcritical  response  with  a  turning  point,  such  as  the  one  shown  in 
Figure  2.6,  nonlinear  structural  parameters  were  chosen  for  the  NACA  64A006  airfoil. 
Care  was  taken  to  choose  parameters  such  that  pitch  and  plunge  excursions  were  still 
valid  in  the  context  of  an  Euler  code,  i.  e.  small  amplitude  LCOs.  With  the  parameters 
chosen,  the  flutter  point,  turning  point,  and  bifurcation  diagram  were  estimated  from 
numerous  computational  runs.  To  demonstrate  the  utility  of  the  stochastic  algorithm  in 
estimating  PDFs  and  probabilities  of  failure  in  this  CFD  application,  a  single  reduced 
velocity  below  the  flutter  point  was  examined  (Vr  =  6.5).  Probability  density  functions  of 
the  LCO  amplitude  in  both  pitch,  aLco,  and  plunge,  x/lco ,  and  of  the  coupled  frequency 
of  the  response,  ujlco ,  estimated  from  the  stochastic  projection  method  via  B-splines  are 
presented.  Refined  B-splinc  surfaces  are  also  presented  to  demonstrate  convergence  of 
the  PDFs  and  thus  provide  a  validation  of  the  stochastic  algorithm. 

5.2  Bifurcation  Diagram  of  the  Aeroelastic  System 

Development  of  the  bifurcation  diagram  begins  with  an  estimation  of  the  flutter 
point.  Recall  from  Sec.  2.4  that  the  flutter  point  is  a  loss  of  linear  stability  and  its  location 
is  unaffected  by  the  nonlinear  parameters.  Gordnier  and  Melville  [7]  efficiently  estimate 
the  flutter  point  of  an  aeroelastic  system  by  first  computing  a  few  cycles  at  reduced 
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(a)  Vr  =  6.7  (b)  Vr  =  6.9 

Figure  5.1  Examples  of  the  pitch  response 

velocities  above  and  below  the  flutter  point  and  then  comparing  amplitude  factors  as  the 
cycles  either  increase  or  decrease.  The  amplitude  factor  is  defined  as  the  ratio  of  the  peak 
magnitude  with  the  magnitude  of  the  previous  peak  with  the  same  sign  (i.e. ,  one  period) 
[7].  The  parameters  given  in  Eq.  (4.40)  with  /3a  =  =  0  provide  a  supercritical  response 

from  which  the  flutter  point  can  be  estimated.  With  an  initial  pitch  of  a(0)  =  0.1°,  cases 
were  run  at  reduced  velocities  of  Vr  —  6.6,  6.7,  6.8,  6.9,  and  7.0.  Examples  of  the 
different  pitch  responses  are  shown  in  Figure  5.1.  Note  that  the  amplitudes  are  of  such 
small  magnitudes  that  any  nonlinear  restoring  forces  would  be  completely  unimportant, 
i.e.  a3,  a5  <C  0.  An  average  of  the  last  few  amplihcation  factors  for  each  of  the  responses 
results  in  the  plot  show  in  Figure  5.2.  From  Figure  5.2  the  flutter  point  is  estimated  to 
occur  at  Vr  ~  6.78. 

Once  the  flutter  point  was  determined,  nonlinear  parameters  were  selected  to  pro¬ 
duce  a  subcritical  response  with  a  turning  point.  Following  the  model  problem,  a  cubic 
spring  constant  of  j3a  <  0  was  chosen  to  destabilize  the  aeroelastic  system  below  the  flut¬ 
ter  point.  The  goal  was  to  choose  a  (3a  such  that  the  aeroelastic  system  was  destabilized 
at  a  small  initial  pitch  angle.  Numerical  experimentation  led  to  a  value  of  f3a  =  —30, 
which  caused  the  aeroelastic  system  to  become  unstable  with  a(0)  ~  3°  at  a  reduced 
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Figure  5.2  Flutter  point  estimation 


Figure  5.3  LCO  due  to  a  subcritical  response,  Vr  =  6.5,  a(0)  =  3°,  j3a  =  —30,  =  500 

velocity  of  Vr  —  6.5.  The  quintic  spring  constant  was  then  chosen  to  restabilizc  the 
system.  Again  through  numerical  experimentation,  a  value  of  7Q  =  500  was  selected. 
This  value  resulted  in  an  LCO  at  Vr  =  6.5  with  an  initial  pitch  angle  of  a(0)  =  3°,  as 
shown  in  Figure  5.3.  For  clarity,  the  parameters  of  the  aeroelastic  system  are  presented 
again.  The  structural  parameters,  with  the  inclusion  of  the  nonlinear  spring  constants, 
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Figure  5.4  Identification  of  the  turning  point 


are 


0.375 

xQ 

=  -0.25 

r2 

a 

=  0.25 

125 

Cv 

=  0.5 

(a 

=  0.5 

0.2 

Pa 

=  -30 

la 

=  500  , 

and  the  aerodynamic  parameters  are 

Poo  =  37612  N/m2  =  0.5489  kg/m3  =  0.85 


(5.1) 


(5.2) 


The  freestream  pressure  and  density  correspond  to  an  altitude  of  25,000  ft. 

With  the  parameters  chosen,  numerous  computer  runs  were  conducted  on  a  range  of 
reduced  velocities  at  and  below  the  flutter  point  to  map  out  the  subcritical  bifurcation. 
In  each  case  a  sufficiently  large  initial  pitch  angle  of  a(0)  =  7°  was  chosen  to  ensure 
the  aeroelastic  system  would  exhibit  LCO.  It  was  found  that  the  system  would  rapidly 
exhibit  LCO  or  rapidly  proceed  to  a  stationary  response  near  the  turning  point,  as  shown 
in  Figure  5.4.  The  turning  point  was  estimated  at  Vr  ~  6.11.  The  resulting  bifurcation 
diagram  is  shown  in  Figure  5.5. 
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Figure  5.5  Bifurcation  diagram  for  the  nonlinear  aeroelastic  system 
5.3  Uncertainty  Quantification 

As  with  the  model  problem,  an  initial  condition  uncertainty  and  a  parametric 
uncertainty  were  added  to  the  formulation.  A  Gaussian  distribution  in  the  uncertainties 
was  assumed.  These  uncertainties  are,  as  before,  given  by 


a(0)  =  a(0)  +fid(0)  (5.3a) 

A*  =  A*  +  6A*  (5.3b) 

For  this  aeroelastic  system,  the  following  mean  and  standard  deviation  values  were  used 

a(0)  =  0  fia  =  -30  (5.4a) 

d(0)  =  1.5°  /3a  =  3  (5.4b) 


Samples  of  the  responses  in  pitch,  plunge,  and  frequency  were  obtained  at  the  nodes 
given  in  Table  3.1  at  a  reduced  velocity  of  Vr  =  6.5.  The  following  is  a  summary  of  the 
results. 
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MCS 


Figure  5.6  B-spline  surfaces  of  the  pitch  response  (/  =  2,  J  =  8) 


I  =  2,  J  =  2,  4,  8 

The  effect  of  increasing  resolution  on  the  £2  axis  (uncertainty  in  the  cubic  spring 
constant)  as  the  resolution  on  the  £1  axis  (uncertainty  in  the  initial  pitch  angle)  remained 
fixed  was  investigated  first.  Each  simulation  was  carried  out  to  tmax  ~  500  or  until  the 
LCO  became  fully  developed.  A  fully  developed  LCO  was  considered  achieved  when  two 
successive  amplitudes,  a^co,  were  within  0.1°  of  each  other.  The  samples  obtained  from 
the  simulations  were  used  to  build  the  B-spline  (piecewise  linear)  surfaces  from  Eq.  (3.3). 
An  example  of  the  resulting  surfaces  is  shown  in  Figure  5.6. 
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Immediately  of  note  from  the  B-spline  surface  in  Figure  5.6  is  that,  as  in  the  model 
problem,  the  symmetric  airfoil  and  the  zero  mean  initial  pitch  angle  lead  to  symmetry 
on  the  £1  axis.  Also,  if  an  LCO  develops  for  a  given  set  of  structural  parameters,  the 
amplitude  of  the  LCO  is  the  same  regardless  of  the  initial  pitch  angle.  These  two  proper¬ 
ties  of  the  B-spline  surface  for  this  aeroelastic  system  can  lead  to  a  reduction  in  both  the 
number  and  the  computational  effort  of  simulations  as  the  resolution  increases  in  either 
axis.  It  will  be  noted  later  in  the  discussion  when  these  properties  were  invoked. 

A  bifurcation  is  observed  on  the  £2  axis.  This  bifurcation  occurs  between  £2  = 
1.53415  and  £2  =  2.5.  This  corresponds  to  /3a  =  —25.398  and  (3a  =  —22.5,  respectively. 
Thus,  as  the  overall  pitch  spring  stiffness  becomes  greater,  i.e.  less  destabilization  from 
the  cubic  spring  constant,  LCO  will  not  develop  for  any  initial  pitch  angle.  This  is  an 
example  of  how  uncertainty  in  the  structural  components  has  an  effect  on  the  LCO  of  an 
aeroelastic  system. 

Also  shown  on  Figure  5.6  are  the  results  of  a  10,000  MCS  obtained  from  the  B-spline 
surface.  (Note:  ten  thousand  randomly  generated  Gaussian  variables  can  occasionally, 
although  rarely,  exceed  the  interval  [—4, 4],  When  this  occurred,  the  value  of  that  random 
variable  was  mapped  to  the  closest  boundary.)  As  expected,  the  MCS  samples  are  con¬ 
centrated  near  the  zero  means.  Since  the  discontinuities  associated  with  the  bifurcations 
are  not  sharply  resolved,  many  realizations  fall  on  values  intermediate  to  the  stationary 
response  and  the  LCO  response.  The  effect  of  these  intermediate  values  can  be  seen  in 
the  PDFs  of  the  response  shown  in  Figure  5.7.  Since  the  bulk  of  intermediate  values  are 
a  result  of  the  poor  resolution  in  the  axis,  increasing  resolution  in  the  £2  axis  had  little 
effect  on  the  the  convergence  of  the  PDF  of  the  pitch  response.  A  very  close  examination 
reveals  that  with  J  =  4  and  J  =  8  the  PDFs  are  essentially  identical.  Further  investi¬ 
gation  of  the  uncertainty  quantification  is,  therefore,  continued  holding  the  resolution  in 
the  £2  axis  fixed  at  J  =  4. 
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(a)  PDFs  of  the  pitch  response 


Figure  5.7 


(b)  Close  up  of  PDFs  of  the  pitch  re¬ 
sponse 

Pitch  Response 


I  =  2,  4,  8;  J  =  4 

Having  examined  the  stochastic  response  surface  with  increasing  resolution  along 
the  £2  axis,  simulations  were  then  performed  with  increasing  resolution  on  the  £1  axis. 
Since  the  amplitudes  of  the  LCOs  were  already  known  from  the  previous  set  of  sim¬ 
ulations,  these  simulations  were  carried  out  just  far  enough  in  time  (imax  <  250)  to 
determine  if  an  LCO  was  developing  or  if  the  response  was  damped.  Examples  of  the 
resulting  B-splinc  surfaces  are  shown  in  Figure  5.8.  Note  that  as  the  resolution  increases, 
the  discontinuity  is  more  sharply  refined.  This  refinement  leads  to  fewer  realizations 
intermediate  to  the  LCO  and  the  stationary  responses. 

As  expected,  fewer  intermediate  realizations  leads  to  a  better  convergence  for  the 
PDFs  of  the  response,  as  shown  in  Figure  5.9.  With  each  increase  in  resolution,  the  PDFs 
show  fewer  intermediate  realizations  while  at  the  same  time  an  area  near  a^co  ~  10° 
remains  unchanged.  Note  that  from  Figure  5.5,  the  mean  responses  at  R  =  6.5  is 
&lco  —  10.1°.  Thus,  from  these  simulations,  it  is  estimated  that  a  bimodal  response 
exists  with  the  primary  response  at  «lco  =  0  and  the  secondary  response  at  «lco  —  10°. 
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Figure  5.8  B-spline  surfaces  of  the  pitch  response  (ctlco) 


(a)  PDFs  of  the  pitch  response 


Figure  5.9 


(b)  Close  up  of  PDFs  of  the  pitch  re¬ 
sponse 

Pitch  Response 


Additionally,  the  uncertainty  bounds  for  the  secondary  response  are  estimated  to  be  from 
otLco  =  9°  to  «lco  =  12°.  The  subcritical  response  and  the  uncertainty  bounds  were 
identified  with  a  sample  size  of  180  simulations. 
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5-4  Refined  B-spline  surface 

A  better  estimate  of  the  subcritical  response  can  easily  be  obtained  by  making  use 
of  the  information  already  gained.  The  locations  of  the  bifurcations  have  been  bounded 
and  can  be  better  refined  with  a  simple  binary  search.  Further  validation  of  this  method 
can  also  be  shown  by  refining  the  B-spline  surface  at  a  different  set  of  nodes.  Thus  the 
nodes  chosen  along  the  £2  axis  were  uniformly  spaced  until  reaching  the  bifurcation  while 
the  nodes  along  the  £1  were  selected  to  refine  the  bifurcations  on  that  axis.  Additionally, 
only  the  positive  half  of  the  £1  axis  was  investigated  and  symmetry  was  assumed  on  the 
negative  half  of  the  axis.  A  few  simulations  were  conducted  to  verify  this  symmetry, 
such  as  the  ones  shown  in  Figure  5.10.  With  the  above  assumptions,  approximately  one 
hundred  simulations  were  performed  to  obtain  the  nodes  for  the  refined  B-spline  shown  in 
Table  5.1.  Since  the  nodes  no  longer  compose  a  regular  mesh,  the  multivariate  B-splinc 
approach  is  difficult  to  implement.  The  10,000  run  MCS  is  accomplished  by  a  simple 
linear  interpolation  scheme,  which  is  equivalent  to  a  piecewise  linear  approximation  of 
the  second  order  B-spline  [30].  A  comparison  of  the  B-spline  surface  obtained  with  1  =  8 
and  J  =  4  with  the  refined  B-spline  surface  and  the  associated  MCS  realizations  is  shown 
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MCS 


(a)  B-spline  surface  from  the  Gaussian  (b)  Refined  B-spline  surface 

nodes  (1  =  8,  J  =  4) 


Figure  5.11  Pitch  Response 


Table  5.1  Nodes  for  the  refined  B-spline  surface 


6 

6 

OlLCO 

6 

6 

OtLCO 

6 

6 

OLLCO 

-4 

-4 

15.6° 

-4 

-1 

10.9° 

-4 

2.03333 

8.89° 

-1.7 

15.6° 

-1.83667 

10.9° 

-3.26 

8.89° 

-1.69667 

0 

-1.83333 

0 

-3.25667 

0 

1.69667 

0 

1.83333 

0 

3.25667 

0 

1.7 

15.6° 

1.83667 

10.9° 

3.26 

8.89° 

4 

15.6° 

4 

10.9° 

4 

8.89° 

-4 

-3 

13.9° 

-4 

0 

10.1° 

-4 

2.06666 

0 

-1.73 

13.9° 

-1.91667 

10.1° 

-3.26 

0 

-1.72667 

0 

-1.91333 

0 

-3.25667 

0 

1.72667 

0 
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(a)  PDFs  of  the  pitch  response 


(b)  Close  up  of  PDFs  of  the  pitch  re¬ 
sponse 


Figure  5.12  Convergence  of  the  pitch  subcritical  response 


in  Figure  5.11.  A  PDF  of  the  pitch  response  obtained  from  the  refined  B-spline  surface 
is  seen  to  achieve  convergence  in  Figure  5.12,  where  the  intermediate  values  between  the 
stationary  response  and  the  LCO  response  are  practically  non-existent. 


Plunge  and  Coupled  Frequency 

The  same  uncertainty  analysis  was  investigated  for  the  plunge  LCO  amplitude  and 
the  coupled  frequency.  B-spline  surfaces  of  the  plunge  LCO  amplitude,  Vlcoi  and  the 
coupled  frequency,  oolco  were  developed  for  the  1  =  8  and  J  =  4  Gaussian  nodes  and 
for  the  refined  nodes.  These  surfaces  are  shown  in  Figures  5.13  and  5.14.  Convergence 
of  the  PDFs  for  these  responses  is  shown  in  Figures  5.15  and  5.16. 


5.5  Probability  of  Failure 

The  PDFs  of  the  LCO  responses  allow  an  estimate  the  probability  of  failure  for  the 
aeroelastic  system.  As  with  the  model  problem,  criteria  need  to  be  chosen  that  define 
the  failure  mode.  Recalling  that  the  subcritical  LCO  response  develops  below  the  flutter 
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(a)  B-spline  surface  from  the  Gaussian  (b)  Refined  B-spline  surface 

nodes  (1  =  8,  J  =  4) 


Figure  5.13  Plunge  Response 


(a)  B-spline  surface  from  the  Gaussian  (b)  Refined  B-spline  surface 

nodes  (1  =  8,  J  =  4) 


Figure  5.14  Frequency  Response 
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(a)  PDFs  of  the  plunge  response 


(b)  Close  up  of  PDFs  of  the  plunge  re¬ 
sponse 


Figure  5.15 


Convergence  of  the  plunge  subcritical  response 
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Figure  5.16  Convergence  of  the  frequency  subcritical  response 
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point,  failure  is  defined  as  the  onset  of  an  LCO.  In  order  to  perform  the  integration  of 
the  PDFs,  and  to  avoid  the  smoothing  error  inherent  in  the  Parzen  windowing  approach, 
failure  in  each  mode  (pitch,  plunge,  and  frequency)  were  set  at 


_  -1  o 

^fail  1 

(5.5a) 

l/fail  =  0.01 

(5.5b) 

^fail  =  1  Hz 

(5.5c) 

With  these  failure  modes  defined,  the  probability  of  failure  was  estimated  from  the  PDFs 
obtained  from  the  Gaussian  nodes  as  well  as  from  the  refined  nodes.  The  results  are 
summarized  in  Table  5.2.  As  resolution  increases,  the  estimation  of  the  probability  of 
failure  improves  until,  at  the  refined  surface,  all  three  probability  of  failures  converge 
to  the  same  value.  The  convergence  of  all  three  probabilities  of  failure  to  the  same 
value  is  an  expected  result  since  the  three  events  (pitch,  plunge,  and  frequency)  are  not 
independent  of  each  other.  If  the  aeroclastic  system  experiences  a  pitch  (plunge)  LCO, 
it  must  also  exhibit  a  plunge  (pitch)  LCO,  and  these  LCOs  will  occur  at  some  frequency. 
Due  to  the  numerous  realizations  at  the  intermediate  values,  the  Gaussian  nodes  provide 
a  conservative  estimate  of  the  probability  of  failure  while  the  refined  nodes  provides 
an  estimate  comparable  to  a  full-order  MCS.  The  total  number  of  actual  simulations 
performed,  however,  are  two  orders  of  magnitude  less  than  that  of  a  10,000  run  MCS. 


Table  5.2  Estimates  of  the  Probability  of  Failure 


Probability  of  Failure  (%) 

Nodes 

Pitch 

Plunge 

Frequency 

7  =  2,  J  =  4 

38.1 

39.7 

36.9 

7  =  4,  J  =  4 

19.6 

20.3 

18.9 

7  =  8,  J  =  4 

10.3 

10.6 

10.1 

Refined 

4.36 

4.36 

4.36 

5-15 


5. 6  Summary 

This  chapter  demonstrated  the  efficiency  of  the  stochastic  projection  method  via 
B-splines  in  estimating  the  probability  of  failure  for  a  CFD  application  of  an  aeroelastic 
system.  The  structure  was  a  low  fidelity,  two  degree  of  freedom,  nonlinear  airfoil  modeled 
in  a  higher  fidelity  inviscid  transonic  flow.  The  number  of  simulations  actually  performed 
was  two  orders  of  magnitude  less  than  a  full  order  MCS.  Verification  of  the  stochastic 
algorithm  was  demonstrated  by  the  convergence  of  the  PDFs  to  a  subcritical  response. 
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VI.  Summary  and  Conclusions 


6. 1  Summary 

A  stochastic  algorithm  based  upon  the  PCE  and  FCE  stochastic  projection  meth¬ 
ods  was  developed  by  employing  multivariate  B-splines.  This  method  is  a  nonintrusive 
(sampling  approach)  that  does  not  require  the  computation  or  storage  of  expected  values 
to  estimate  the  stochastic  projection  of  the  response.  Instead,  the  B-spline  approach  is 
a  collocation  method  that  fits  a  surface  to  the  nodes  provided.  The  nodes  were  selected 
both  by  a  Gaussian  distribution  as  well  as  a  uniform  distribution  of  probabilities. 

An  aeroelastic  model  problem  was  first  investigated  by  a  traditional  Monte  Carlo 
approach.  Bifurcation  diagrams,  PDFs,  and  probabilities  of  failure  were  provided  for 
comparison  with  the  stochastic  algorithm.  The  computational  effort  for  this  model  prob¬ 
lem  was  shown  to  be  quite  large. 

A  small  sample  set  based  upon  a  Gaussian  distribution  of  probabilities  was  then 
obtained.  A  B-splinc  surface  was  developed  from  which  a  fast  MGS  could  be  performed. 
Although  the  LCOs  developed  were  due  to  initial  conditions  and  large  time  integrations, 
the  sampling  approach  removes  the  samples  from  the  time  domain  by  tracking  the  LCO 
amplitude  once  the  LCO  is  fully  devloped.  Results  from  the  stochastic  algorithm  were 
in  excellent  agreement  with  the  traditional  MGS  approach.  Additionally,  the  stochastic 
algorithm  was  two  orders  of  magnitude  more  efficient  than  the  MGS  approach. 

The  applicability  of  the  stochastic  algorithm  was  investigated  for  the  more  compu¬ 
tationally  demanding  CFD  simulation  of  an  aeroelastic  system.  An  inviscid  aeroelastic 
code  (ELILER-AE)  was  developed,  verified,  and  validated  to  perform  this  investigation. 
Parameters  were  selected  to  obtain  a  subcritical  response  with  a  turning  point.  This 
response  is  considered  high  risk  since  large  amplitude  LCOs  can  develop  below  the  clas¬ 
sically  defined  flutter  point.  Monte  Carlo  simulations  were  impractical  for  this  problem 
since  each  simulation  took  36  to  48  hours  to  perform.  Validity  of  the  stochastic  projection 
method  via  B-splines  was  determined  by  convergence  of  the  PDFs  and  the  probability 
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of  failure.  Again,  efficiency  was  two  orders  of  magnitude  greater  than  that  of  the  Monte 
Carlo  approach,  making  this  stochastic  algorithm  applicable  to  large  time,  CFD  applica¬ 
tions. 

6.2  Conclusions 

Based  upon  this  research  effort,  the  following  conclusions  were  reached: 

1.  For  large  time  and  highly  nonlinear  problems,  nonintrusive  (sampling)  methods  are 
more  efficient  than  intrusive  methods.  This  point  was  alluded  to  in  the  introduction. 
For  large  time  problems,  the  frequency  content  in  the  stochastic  domain  at  any 
particular  time  can  become  so  large,  or  contain  so  many  discontinuities,  that  even 
a  basis  that  might  be  considered  optimal  would  require  such  a  high  order  expansion 
that  all  computational  efficiency  is  lost. 

2.  Bifurcations  are  represented  by  discontinuities  in  the  stochastic  domain,  making 
spectral  approaches  inefficient.  When  the  spectral  approach  is  recognized  as  an 
attempt  to  approximate  a  global  response  surface,  its  limitations  are  the  same  as 
in  any  interpolation  problem.  An  infinite  expansion  is  required  to  sharply  define  a 
discontinuity. 

3.  Efficiency  can  be  regained  by  distribution  based  sampling  and  approximating  the 
global  response  surface  with  a  multivariate  B-spline.  Splines  are  the  most  appropri¬ 
ate  method  for  approximating  a  peicewise  continuous  function  in  any  interpolation 
problem.  The  selection  of  nodes  should  initially  be  based  on  the  distribution  of 
input  uncertainty.  Refinement  comes  once  information  is  obtained  from  an  initial 
set  of  simulations. 

4.  The  stochastic  projection  method  via  B-splines  is  applicable  to  CFD  based  research. 
Uncertainty  quantification  is  an  expensive  proposition.  As  CFD  codes  become 
faster  due  to  improvements  of  algorithms  and  processor  power,  the  ability  to  con¬ 
duct  multiple  simulations,  though  not  at  the  level  of  a  full-order  MCS,  moves  into 


6-2 


the  realm  of  possibilities.  The  stochastic  algorithm  presented  here  makes  uncer¬ 
tainty  quantification  possible  at  the  current  level  of  processor  power. 

6.3  Recommendations  for  Future  Research 

The  stochastic  algorithm  was  applied  to  an  inviscid,  transonic  aeroelastic  prob¬ 
lem.  The  next  step  in  the  study  of  aeroelastic  responses  would  be  to  include  viscous 
fluxes,  turbulence  modeling,  and  eventually  a  three  dimensional  configuration.  Since  the 
method  is  nonintrusive,  any  research  or  commercial  code  could  be  used  to  conduct  this 
investigation. 

Future  research  with  the  stochastic  algorithm  should  include: 

1.  The  inclusion  of  higher  order  multivariate  splines  The  piecewise  linear  interpolation 
was  nearly  optimal  for  the  two  dimensional  aeroelastic  system.  This  will  not  always 
be  the  case.  Indeed,  for  the  projection  with  respect  to  £2,  Figures  3.2(b)  and  5.11(b), 
a  higher  order  interpolation  would  provide  a  better  fit  than  the  piecewise  linear 
interpolation.  However,  the  bifurcation  in  Figure  5.11(b)  produces  overshoots  with 
higher  order  B-splincs.  The  Akima  spline  [30]  is  shape  preserving  and  may  be  more 
appropriate  for  this  type  of  projection. 

2.  Determination  of  the  maximum  number  of  uncertainties  that  can  be  included  in  the 
algorithm.  The  number  of  simulations  increases  geometrically  if  the  same  number 
of  nodes  are  used  in  each  dimension  of  uncertainty.  Thus,  there  is  a  practical  limit 
to  the  number  of  uncertainties  that  can  be  modeled  before  the  computational  effort 
reaches  that  of  a  Monte  Carlo  simulation.  Investigating  each  dimension  individually 
and  determining  the  maximum  nodes  required  before  moving  to  the  next  dimension, 
as  was  done  with  the  CFD  application,  can  lead  to  an  algebraic  increase  in  effort 
as  dimensions  are  added. 

3.  Multiple  distributions  of  uncertainties.  In  the  CFD  application,  both  Gaussian 
nodes  and  uniformly  spaced  nodes  were  used  to  build  the  B-spline  surface.  This 
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is  an  indication  that,  unlike  other  stochastic  projection  methods,  the  distribution 
of  uncertainties  do  not  need  to  be  homogeneous  among  the  random  variables.  The 
intervals  and  nodes  can  be  selected  based  upon  the  individual  distribution  of  the 
random  variables. 

4.  Automation.  Once  the  Gaussian  nodes  were  computed,  an  interactive  approach 
was  employed  to  refine  the  sampling  mesh.  Automating  this  process  would  be  a 
major  undertaking  but  could  provide  a  great  value  to  the  aircraft  designer. 

5.  Dimensional  Analysis  It  was  demonstrated  that  different  LCO  solutions  were  ob¬ 
tained  from  the  matched  and  non-matched  point  assumptions.  The  proper  nondi- 
mensional  parameters  for  the  fluid-structure  coupling  should  be  determined. 
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Appendix  A.  Coupling  of  the  Nondimensional  Parameters 

As  early  as  1973  [50]  and  as  recently  as  2004  [51]  it  has  been  recognized  that  if  Mach 
number  is  specified  when  solving  an  aeroelastic  system,  compatibility  should  be  enforced 
between  freestream  Mach  number,  freestream  density,  and  the  reduced  velocity.  This 
compatibility  is  referred  to  as  a  matched  point  analysis.  Other  researchers  [34,  40]  prefer 
a  non-matched  point  solution  where  only  the  nondimensional  parameters  are  specified 
and  it  is  assumed  that  a  similarity  solution  which  conforms  to  a  matched  point  solution 
exists.  It  will  be  shown  that  both  approaches  have  limitations  that  exist  due  to  the 
coupling  of  the  fluid  and  structure. 

A  straightforward  dimensional  analysis  of  the  Euler  equations  leads  to  the  nondi¬ 
mensional  parameters  shown  in  Eq.  (4.1).  Of  particular  note  is  the  nondimensionalization 
of  the  farfield  boundary  condition  (here  the  asterisk  indicates  a  nondimensional  quantity) 


*  1 

Poo  =  1 

(A.  la) 

ulc  =  cos  a 

(A. lb) 

V*oo  =  sin« 

(A.lc) 

Po°  ~  7  Ml 

(A.  Id) 

Thus,  a  properly  nondimensionalized  Euler  code  only  requires  the  freestream  Mach  num¬ 
ber  and  angle  of  attack  to  compute,  say,  a  pressure  distribution  over  an  airfoil.  The 
resulting  solution  is  a  similarity  solution  for  any  two  given  thermodynamic  quantities, 
say  pressure  (poo)  and  density  (poo). 

A  similar  dimensional  analysis  can  be  performed  on  the  pitch  and  plunge  airfoil 
equations  with  modeled  aerodynamics,  Eq.  (2.1).  Lift  and  moment  are  scaled  by  the 
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dynamic  pressure  and  chord,  that  is 


r=y*coV£cC,(t)  (A. 2a) 

M  =  £c2Cm(t)  (A. 2b) 

It  is  important  to  recognize  that  for  the  modeled  aerodynamics,  the  lift  and  moment 
coefficients  are  functions  only  of  time  and  the  flow  is  assumed  incompressible.  Thus,  the 
nondimensionalization  that  leads  to  the  reduced  velocity  (Vr),  frequency  ratio  (our),  and 
mass  ratio  (/xs)  provide  a  similarity  solution  for  any  given  dynamic  pressure  ( p ^  and  I4o) 
and  structural  parameters  (m,  ua,  lvh,  xa,  and  ra). 

When  coupling  the  Euler  equations  to  the  structural  equations,  a  coupling  also 
occurs  in  the  nondimensional  parameters.  The  reduced  velocity,  for  example,  is  given  by 
(for  convenience,  all  lengths  in  this  section  are  scaled  by  the  chord  c) 

Vr  =  —  (A.  3) 

u JaC 

Typically  the  reduced  velocity  is  taken  as  a  free  parameter  in  the  structural  problem,  but 
it  will  be  shown  that  for  the  coupled  fluid-structure  problem  the  reduced  velocity  should 
instead  be  derived  from  the  remaining  system  parameters. 

To  demonstrate  this  assertion,  consider  the  definition  of  Mach  number 

=  —  (A. 4) 

^oo 

where  a ^  is  the  freestream  speed  of  sound  and  is  given  by 

(A.5) 
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(a)  Matching  density  with  remaining  (b)  Comparison  with  standard  day 

parameters 


Figure  A.l  Density  determined  by  varying  uja  and  p^ 


The  reduced  velocity  is  then  rewritten  as 


Vr 


Afoo  I  Po c 

^ a C  V  P°° 


(A. 6) 


Rearranging  terms  to  solve  for  freestream  density  yields 


ipoqML 

c o2ac2v;2 


(A.  7) 


Equation  (A. 7)  shows  that  for  a  given  Mach  number  and  reduced  velocity  the  freestream 
density  is  a  function  of  freestream  pressure  (a  fluid  parameter)  and  of  the  uncoupled 
natural  frequency  in  pitch  (a  structural  parameter).  The  implications  of  this  relationship 
between  fluid  and  structure  parameters  is  shown  in  Figure  A.l.  To  obtain  these  results, 
the  Mach  number  was  selected  as  M ^  =  0.85  and  the  reduced  velocity  as  Vr  =  11. 
The  freestream  pressure  was  allowed  to  vary  with  altitude  according  to  standard  day 
tables  and  the  uncoupled  natural  frequency  in  pitch  was  allowed  to  vary  from  ua  =  10 
to  tua  =  50.  Figure  A. 1(a)  shows  the  large  variation  obtained  in  density,  particularly 
near  sea  level  pressure.  These  results  are  typical  of  a  non-matched  point  analysis  where 
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only  the  nondimensional  parameters  (to  include  Mach  number  and  reduced  velocity)  are 
specified.  In  order  to  obtain  the  solution  for  ojq  =  10 Hz  at  sea  level  pressure,  density  is 
on  the  order  of  8  kg/m3,  an  unrealistic  value  for  air.  Even  at  low  pressures,  p^  ~  20000 
N/m2  (approximately  40,000  ft  pressure  altitude),  density  varies  from  greater  than  sea 
level  values  (poo  ~  1-67  kg/m3)  to  approximately  70,000  ft  density  altitude.  In  contrast, 
a  matched  point  solution  specifies  a  pressure  and  density  that  vary  with  a  standard  day. 
Figure  A. 1(b)  shows  where  the  standard  day  pressure  and  density  lie  with  comparison  to 
the  non-matched  point  data.  Note  that  the  matched  point  solution  allows  only  a  small 
range  of  possible  structural  parameters.  One  point  is  clearly  evident  from  Figures  A. 1(a) 
and  (b),  i.e. ,  the  matched  point  and  non-matched  point  parameters  will  not  necessarily 
provide  the  same  solution. 

The  above  analysis  indicates  similarity  solutions  for  a  general  structure  and  general 
freestream  conditions  are  difficult  to  obtain.  In  order  to  obtain  solutions  of  realistic 
structures  in  realistic  freestream  conditions,  the  following  parameters  should  be  specified: 
two  thermodynamic  properties  (such  as  p^  and  p^)  and  the  structural  parameters  (m, 
u>a,  iOh,  xQ  and  ra).  The  reduced  velocity  and  frequency  and  mass  ratios  become  derived 
quantities.  The  only  free  parameter  is  then  Mach  number. 
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Appendix  B.  A  Summary  of  Computer  Runs 

The  following  are  the  contents  of  the  data  compact  disc  (CD). 

1.  Directory  containing  EULER- AE  source  files  and  executable.  Also  contains  a  make¬ 
file  for  a  64-bit  processor.  Two  grid  files  are  included,  one  for  the  NACA  64A006 
and  one  for  NACA  64A010  (asymmetric).  The  input  file  and  the  script  to  submit 
the  job  are  also  included. 

2.  Directory  containing  the  time  histories  of  the  cases  run  to  generate  the  results  in 
this  dissertation. 

3.  Directory  containing  restart  files  for  some  of  the  cases  that  were  run. 

A  brief  explanation  of  how  to  use  the  ELILER-AE  code  is  presented  along  with  a 
summary  of  the  computer  runs. 

The  Linux/UNIX  version  of  EULER- AE  is  designed  to  be  run  from  the  command 
line  with  the  following  statement: 

,/euler_ae  < input _file_name>  <grid_f ile_name>  <restart_f ile_name> 

A  sample  of  the  input  file  is  shown  in  Figure  B.l.  The  Start/Restart  information 
provides  names  for  the  output  files.  The  input  files  restartfile  and  gridfile  are 
ignored  since  these  are  obtained  from  the  command  line.  The  name  of  the  restart  file  to 
be  read  in  subsequent  runs  is  saved  as  binary  in  newrestartf  ile.  The  final  flowfield, 
at  the  end  of  steady  state  convergence  or  the  last  time  accurate  solution,  is  saved  as 
ASCII  and  Tecplot™  format  in  fflow.  The  animation  file  for  Tecplot™  is  saved  as 
ASCII  in  inflow.  The  time  history  file,  timehist,  is  saved  in  ASCII.  The  first  variable  is 
nondimensional  time,  the  second  is  pitch,  the  third  is  plunge,  the  fourth  is  the  coefficient 
of  lift,  and  the  last  is  the  moment  coefficient.  The  grid  size  imax  and  jmax  must  match 
the  grid  input  file. 

The  inflow  and  reference  length  parameters  are  entered  next.  The  code  expects 
dimensional  SI  units  for  input  and  will  perform  the  nondimensionalizations.  Some  of 
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Input  paramters  for  EULER_AE 


St art /Re st art  Information 

1  !  (restart)  0  =  restart,  1  =  new  case 
restartcase293.dat  !  (restartf ile)  Name  of  restart  file  to  be  read 
restartcase293.dat  !  (newrestartf ile)  Name  of  new  restart  file  to  be  written 
fflowcase293.dat  !  (fflow)  Final  flow  file 
mflowcase293.dat  !  (mflow)  Movie  flow  file  (animation) 
timehistcase293.dat  !  (timehist)  Time  history  of  pitch,  plunge,  C_l,  and  C_m 
naca64a006v3.dat  !  (gridfile)  Name  of  grid  file 
257  !  (imax) 

65  !  (jmax) 

Inflow  and  reference  length 

37612.0  !  (p_in)  Free-stream  pressure  (N/nT2) 

0.5489  !  (rho_in)  Free-stream  density  (kg/nT3) 

1.7887e-5  !  (mu_in)  Free-stream  viscosity  (kg/ms''2) 

287.05  !  (R_gas)  Gas  constant  (N*m/kg*K) 

1.4  !  (gamma)  Ratio  of  specific  heats 
0.850  !  (M_in)  Free-stream  Mach 
6.00  !  (aoa_deg)  Angle  of  attack  (degrees) 

0.72  !  (Pr)  Prandtl  Number 

1.0000  !  (ref_len)  Reference  length  (m) 

288.15  !  (T_w)  Wall  Temperature  (K)  (Make  this  negative  for  adiabatic  wall) 

Program  Parameters 

25000  !  (mmax)  maximum  number  of  time  steps  (steady-state) 

26000  !  (jump)  output  interval 
1.0e-5  !  (tol)  Convergence  tolerance  (steady  state) 

500.0  !  (tend)  Final  time  (time  accurate) 

1  !  (istep)  0  =  use  computed  time  step,  1  =  user  defined  time  step  (time  accurate) 

0.00020  !  (user_dt)  user  defined  time  step  (time  accurate) 

Solution  Parameters 

0  !  (march)  0  =  time  explicit,  1  =  dual  time  implicit 
1.00  !  (CFL)  Stability  time  step  control 
1  !  (method)  0  =  JST  1  =  Roe 
1.0  !  (psi)  0  =  1st  order  Roe,  1  =  higher  order  Roe 
0.333333333333333  !  (kappa)  -1  =  0(2)  upwind,  0  =  0(2)  Fromm,  1/3  =  0(3)  upwind,  1  =  0(2)  cen.  dif f . 
1  !  (limiter)  0  =  no  limiter,  1  =  minmod,  2  =  superbee 
0.0  !  (ent)  0  =  no  entropy  fix,  1  =  Harten  and  Hyman,  4  =  Kermani  and  Piet 

Boundary  Conditions 

0  !  (iminbndry)  imin  boundary  condition,  0  =  farfield,  1  =  wall 

0  !  (imaxbndry)  imax  boundary  condition,  0  =  farfield,  1  =  wall 

1  !  (jminbndry)  jmin  boundary  condition,  0  =  farfield,  1  =  wall 

0  !  (jmaxbndry)  jmax  boundary  condition,  0  =  farfield,  1  =  wall 

65  !  (icutl)  Set  to  0  unless  using  a  C-grid 

193  !  (icut2)  Will  be  set  to  imax  +  1  if  icutl  =  0 

Airfoil  Parameters 

0.375  !  (x_cg)  Distance  (°/0  chord)  from  the  airfoil  leading  edge  to  the  CG 

-0.25  !  (x_alpha)  Distance  (°/«  semi-chord)  from  the  CG  to  the  elastic  axis 

0.25  !  (r_alpha_sq)  Radius  (°/«  semi-chord)  of  gyration  (squared) 

0.2  !  (omega_bar)  Frequency  ratio 
0.5  !  (zeta_h) 

0.5  !  (zeta_alpha) 

-25.398  !  (beta_alpha) 

500 . 0  !  (gamma_alpha) 

125.0  !  (m_bar)  mass  ratio 

6.50  !  (V_bar)  Reduced  velocity  (nondimens ional)  (°/0  semi-chord) 

0.00  !  (plunge)  Initial  plunge  (nondimens ional)  (positive  in  the  up  direction) 


Figure  B.l  Typical  input  file  for  EULER-AE 
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these  parameters  (mu_in,  Pr,  T_w)  are  for  a  viscous  code  and  are  ignored  in  the  current 
version  of  EULER- AE. 

The  program  parameters  control  the  number  of  steady  state  iterations  and  the 
maximum  time  for  the  time  accurate  integration.  For  a  steady  state  solution  only,  set 
ramax  to  the  desired  number  of  iterations  and  tend  to  0.  To  start  a  time  accurate  solution 
from  a  converged  steady  state  solution  or  to  restart  a  time  accurate  solution,  set  restart 
to  1,  mmax  to  0  and  tend  to  the  desired  maximum  time.  Since  the  restart  hie  saves  the 
time  the  last  integration  was  performed  to,  tend  must  be  greater  that  that  time  in  order 
for  the  program  to  continue  integrating.  The  output  interval  jump  controls  how  often 
each  animation  hie  is  written  to  mf  low. 

The  solution  parameters  control  the  type  of  time  integration  to  use  and  the  scheme 
used  to  compute  the  hux  vectors.  The  dual  time  step  option  (march  set  to  1)  is  not  fully 
functional  and  should  not  be  used.  The  code  is  capable  of  computing  the  hux  vectors  from 
the  Jameson-Schmidt-Turkel  scheme  or  the  Roe  scheme.  The  Roe  scheme  incorporates 
the  MUSCL  scheme  to  allow  for  higher  order  spatial  integration.  Additionally,  two 
different  limiters  and  entropy  correction  schemes  can  be  chosen.  The  parameters  shown 
in  Figure  B.l  were  used  throughout  this  study,  i.e.  a  second  order  explicit  time  integration 
with  a  third  order  Roe  scheme  with  the  minmod  limiter  and  no  entropy  correction. 

Two  different  boundary  conditions  are  set  explicitly  -  farheld  and  wall.  The  third 
is  selected  based  upon  the  value  of  icutl.  This  is  the  cut  condition.  If  icutl  is  set  to 
0,  no  cut  boundaries  are  assumed  (e.g.  a  supersonic  ramp).  For  the  257  x  65  point  grid 
used  in  this  study,  the  values  of  icutl  and  icut2  shown  in  Figure  B.l  indicate  a  cut 
boundary  condition  exists  from  1  to  65  and  193  to  257. 

Finally,  the  structural  parameters  of  the  airfoil  are  entered.  These  parameters  are 
xcg  (x_cg),  xa  (x_alpha),  rjj,  (r_alpha_sq),  ur  (omega_bar),  (h  (zeta_h),  (a  (zeta_alpha), 
f3a  (beta_alpha),  7a  (gamma_alpha),  and  /is  (m_bar).  The  reduced  velocity  Vr  (V_bar) 
and  the  initial  nondimensional  plunge  displacement  y(0)  (plunge)  are  also  entered  here. 
Note  the  initial  rates  d(0)  and  y(0)  are  assumed  to  be  zero. 
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The  parameters  that  changed  from  run  to  run  were  «(())  (aoa_deg),  f3a  (beta_alpha), 
7a  (gamma_alpha) ,  and  Vr  (V_bar).  Table  B.l  summarizes  each  of  the  cases.  Cases  1,  2, 
and  3  are  not  shown  since  they  were  used  to  perform  final  verification  of  the  code. 
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Table  B.l  Summary  of  computer  runs 


Case 

a(0) 

/3a 

7  a 

vr 

4 

1.0 

-3.0 

0.0 

6.0 

5 

2.0 

-3.0 

0.0 

6.0 

6 

3.0 

-3.0 

0.0 

6.0 

7 

4.0 

-3.0 

0.0 

6.0 

8 

5.0 

-3.0 

0.0 

6.0 

9 

6.0 

-3.0 

0.0 

6.0 

10 

7.0 

-3.0 

0.0 

6.0 

11 

1.0 

-10.0 

0.0 

6.6 

12 

2.0 

-10.0 

0.0 

6.6 

13 

3.0 

-10.0 

0.0 

6.6 

14 

4.0 

0.0 

6.6 

15 

5.0 

0.0 

6.6 

16 

6.0 

0.0 

6.6 

17 

7.0 

0.0 

6.6 

18 

1.0 

0.0 

6.5 

19 

2.0 

0.0 

6.5 

20 

3.0 

0.0 

6.5 

21 

4.0 

0.0 

6.5 

22 

5.0 

0.0 

6.5 

23 

6.0 

0.0 

6.5 

24 

7.0 

6.5 

25 

1.0 

6.5 

26 

2.0 

6.5 

27 

3.0 

2000.0 

6.5 

28 

4.0 

2000.0 

6.5 

29 

5.0 

6.5 

30 

6.0 

6.5 

31 

7.0 

6.5 

32 

1.0 

6.5 

33 

2.0 

6.5 

34 

3.0 

6.5 

35 

4.0 

6.5 

36 

5.0 

2500.0 

6.5 

37 

6.0 

2500.0 

6.5 

38 

7.0 

2500.0 

6.5 

39 

1.0 

6.5 

40 

2.0 

6.5 

41 

3.0 

2000.0 

6.5 

42 

4.0 

2000.0 

6.5 

43 

5.0 

2000.0 

6.5 
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Table  A.l  Summary  of  computer  runs  (continued) 


Case 

a(0) 

Pa 

la 

Vr 

44 

6.0 

-40.0 

2000.0 

6.5 

45 

7.0 

-40.0 

2000.0 

6.5 

46 

1.0 

-40.0 

1500.0 

6.5 

47 

2.0 

-40.0 

1500.0 

6.5 

48 

3.0 

-40.0 

1500.0 

6.5 

49 

4.0 

-40.0 

1500.0 

6.5 

50 

5.0 

-40.0 

1500.0 

6.5 

51 

6.0 

-40.0 

1500.0 

6.5 

52 

7.0 

-40.0 

1500.0 

6.5 

53 

1.0 

-40.0 

1000.0 

6.5 

54 

2.0 

-40.0 

1000.0 

6.5 

55 

3.0 

-40.0 

1000.0 

6.5 

56 

4.0 

-40.0 

1000.0 

6.5 

57 

5.0 

-40.0 

1000.0 

6.5 

58 

6.0 

-40.0 

1000.0 

6.5 

59 

7.0 

-40.0 

1000.0 

6.5 

60 

1.0 

-30.0 

500.0 

6.5 

61 

2.0 

-30.0 

500.0 

6.5 

62 

3.0 

-30.0 

500.0 

6.5 

63 

4.0 

-30.0 

500.0 

6.5 

64 

5.0 

-30.0 

500.0 

6.5 

65 

6.0 

-30.0 

500.0 

6.5 

66 

7.0 

-30.0 

500.0 

6.5 

67 

7.0 

-30.0 

500.0 

6.0 

68 

7.0 

-30.0 

500.0 

6.1 

69 

7.0 

-30.0 

500.0 

6.2 

70 

7.0 

-30.0 

500.0 

6.3 

71 

7.0 

-30.0 

500.0 

6.4 

72 

7.0 

-30.0 

500.0 

6.6 

73 

7.0 

-30.0 

500.0 

6.12 

74 

7.0 

-30.0 

500.0 

6.14 

75 

7.0 

-30.0 

500.0 

6.16 

76 

7.0 

-30.0 

500.0 

6.18 

77 

7.0 

-30.0 

500.0 

6.105 

78 

7.0 

-30.0 

500.0 

6.11 

79 

7.0 

-30.0 

500.0 

6.115 

80 

2.25 

-30.0 

500.0 

6.5 

81 

2.50 

-30.0 

500.0 

6.5 

82 

2.75 

-30.0 

500.0 

6.5 

83 

6.0 

-27.0 

500.0 

6.5 
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Table  A.l  Summary  of  computer  runs  (continued) 


Case 

a(0) 

Pa 

7a 

Vr 

84 

6.0 

-24.0 

6.5 

85 

6.0 

-21.0 

6.5 

86 

6.0 

-18.0 

500.0 

6.5 

87 

6.0 

-33.0 

500.0 

6.5 

87b1 

6.0 

-33.0 

500.0 

6.5 

88 

6.0 

-36.0 

6.5 

88b 

6.0 

-36.0 

6.5 

89 

6.0 

-39.0 

6.5 

89b 

6.0 

-39.0 

6.5 

90 

6.0 

-42.0 

6.5 

91 

7.0 

6.65 

92 

7.0 

6.7 

93 

7.0 

6.75 

94 

2.8 

6.5 

95 

2.85 

6.5 

96 

2.9 

6.5 

97 

2.95 

6.5 

98 

7.0 

6.55 

99 

7.0 

6.6 

7.0 

6.65 

6.0 

-23.5 

6.5 

6.0 

-23.0 

6.5 

103 

6.0 

-22.5 

500.0 

6.5 

104 

6.0 

-22.0 

500.0 

6.5 

105 

6.0 

-23.9 

500.0 

6.5 

6.0 

-23.8 

6.5 

6.0 

-23.7 

6.5 

6.0 

-23.6 

6.5 

7.0 

30.0 

6.8 

7.0 

30.0 

6.9 

111 

7.0 

30.0 

7.0 

112 

7.0 

30.0 

7.1 

113 

2.86 

-30.0 

6.5 

114 

2.865 

-30.0 

6.5 

115 

2.87 

-30.0 

6.5 

116 

2.875 

-30.0 

6.5 

116b 

2.875 

-30.0 

6.5 

117 

2.88 

-30.0 

6.5 

117b 

2.88 

-30.0 

6.5 

118 

2.885 

-30.0 

6.5 

1  “b” 

indicates  a 

restart  file 
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Table  A.l  Summary  of  computer  runs  (continued) 


Case 

a(0) 

Pa 

7a 

Vr 

119 

2.89 

-30.0 

500.0 

6.5 

120 

-2.875 

-30.0 

500.0 

6.5 

121 

-2.87 

-30.0 

500.0 

6.5 

122 

0.5 

-42.0 

500.0 

6.5 

123 

1.0 

-42.0 

500.0 

6.5 

124 

1.5 

-42.0 

500.0 

6.5 

125 

2.0 

-42.0 

500.0 

6.5 

126 

2.5 

-42.0 

500.0 

6.5 

127 

3.0 

-42.0 

500.0 

6.5 

128 

3.0 

-39.0 

500.0 

6.5 

129 

3.0 

-36.0 

500.0 

6.5 

130 

3.0 

-33.0 

500.0 

6.5 

131 

3.0 

-30.0 

500.0 

6.5 

132 

7.0 

30.0 

500.0 

7.1 

133 

7.0 

30.0 

500.0 

7.2 

134 

2.95 

-42.0 

500.0 

6.5 

135 

2.9 

-42.0 

500.0 

6.5 

136 

2.85 

-42.0 

500.0 

6.5 

137 

2.8 

-42.0 

500.0 

6.5 

138 

2.75 

-42.0 

500.0 

6.5 

139 

2.7 

-42.0 

500.0 

6.5 

140 

2.69 

-42.0 

500.0 

6.5 

141 

2.68 

-42.0 

500.0 

6.5 

142 

2.67 

-42.0 

500.0 

6.5 

143 

2.66 

-42.0 

500.0 

6.5 

144 

2.65 

-42.0 

500.0 

6.5 

145 

2.64 

-42.0 

500.0 

6.5 

146 

2.63 

-42.0 

500.0 

6.5 

147 

2.62 

-42.0 

500.0 

6.5 

148 

2.61 

-42.0 

500.0 

6.5 

149 

2.60 

-42.0 

500.0 

6.5 

150 

2.59 

-42.0 

500.0 

6.5 

151 

7.0 

30.0 

500.0 

7.3 

151b 

7.0 

30.0 

500.0 

7.3 

152 

7.0 

30.0 

500.0 

7.4 

152b 

7.0 

30.0 

500.0 

7.4 

153 

7.0 

30.0 

500.0 

7.5 

153b 

7.0 

30.0 

500.0 

7.5 

154 

7.0 

30.0 

500.0 

7.2 

154b 

7.0 

30.0 

500.0 

7.2 
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Table  A.l  Summary  of  computer  runs  (continued) 


Case 

a(0) 

Pa 

7a 

Vr 

155 

2.585 

-42.0 

500.0 

6.5 

156 

2.58 

-42.0 

500.0 

6.5 

157 

2.575 

-42.0 

500.0 

6.5 

158 

2.57 

-42.0 

500.0 

6.5 

159 

2.565 

-42.0 

500.0 

6.5 

160 

2.56 

-42.0 

500.0 

6.5 

161 

2.555 

-42.0 

500.0 

6.5 

162 

2.55 

-42.0 

500.0 

6.5 

163 

2.545 

-42.0 

500.0 

6.5 

164 

2.54 

-42.0 

500.0 

6.5 

165 

2.535 

-42.0 

500.0 

6.5 

166 

2.53 

-42.0 

500.0 

6.5 

167 

2.525 

-42.0 

500.0 

6.5 

168 

-2.545 

-42.0 

500.0 

6.5 

169 

-2.55 

-42.0 

500.0 

6.5 

170 

2.61 

-39.0 

500.0 

6.5 

171 

2.62 

-39.0 

500.0 

6.5 

172 

2.63 

-39.0 

500.0 

6.5 

173 

2.64 

-39.0 

500.0 

6.5 

174 

2.65 

-39.0 

500.0 

6.5 

175 

2.66 

-39.0 

500.0 

6.5 

176 

2.67 

-39.0 

500.0 

6.5 

177 

2.68 

-39.0 

500.0 

6.5 

178 

2.69 

-39.0 

500.0 

6.5 

179 

2.555 

-39.0 

500.0 

6.5 

180 

2.56 

-39.0 

500.0 

6.5 

181 

2.565 

-39.0 

500.0 

6.5 

182 

2.57 

-39.0 

500.0 

6.5 

183 

2.575 

-39.0 

500.0 

6.5 

184 

2.58 

-39.0 

500.0 

6.5 

185 

2.585 

-39.0 

500.0 

6.5 

186 

2.59 

-39.0 

500.0 

6.5 

187 

2.595 

-39.0 

500.0 

6.5 

188 

2.6 

-39.0 

500.0 

6.5 

189 

2.605 

-39.0 

500.0 

6.5 

190 

2.66 

-36.0 

500.0 

6.5 

191 

2.665 

-36.0 

500.0 

6.5 

192 

0.5 

30.0 

500.0 

6.8 

192b 

0.5 

30.0 

500.0 

6.8 

193 

0.5 

30.0 

500.0 

6.9 
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Table  A.l  Summary  of  computer  runs  (continued) 


Case 

a(0) 

Pa 

7a 

Vr 

194 

7.0 

195 

7.1 

196 

500.0 

6.85 

197 

2.67 

-36.0 

500.0 

6.5 

198 

2.675 

-36.0 

500.0 

6.5 

199 

6.75 

2.76 

-33.0 

6.5 

201 

2.765 

-33.0 

6.5 

202 

2.77 

-33.0 

6.5 

203 

2.775 

-33.0 

6.5 

204 

2.9 

-27.0 

6.5 

205 

2.925 

-27.0 

6.5 

206 

2.95 

-27.0 

500.0 

6.5 

207 

2.975 

-27.0 

500.0 

6.5 

208 

-27.0 

500.0 

6.5 

209 

2.755 

-33.0 

6.5 

2.75 

-33.0 

6.5 

211 

2.745 

-33.0 

6.5 

212 

2.74 

-33.0 

6.5 

213 

2.735 

-33.0 

6.5 

214 

2.73 

-33.0 

6.5 

215 

-27.0 

6.5 

216 

3.10 

-27.0 

500.0 

6.5 

217 

3.15 

-27.0 

500.0 

6.5 

218 

3.20 

-27.0 

500.0 

6.5 

219 

3.25 

6.5 

3.30 

-27.0 

6.5 

221 

3.35 

-27.0 

500.0 

6.5 

222 

■n 

30.0 

500.0 

6.7 

223 

BH 

30.0 

500.0 

6.8 

224 

msm 

6.9 

225 

■  I . 

7.0 

226 

3.055 

-27.0 

6.5 

227 

3.06 

-27.0 

6.5 

228 

3.065 

-27.0 

6.5 

229 

-27.0 

6.5 

3.075 

-27.0 

6.5 

231 

-27.0 

6.5 

232 

3.1 

-23.9 

6.5 

233 

3.125 

-23.9 

6.5 
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Table  A.l  Summary  of  computer  runs  (continued) 


Case 

a(0) 

Pa 

7  a 

Vr 

234 

3.15 

-23.9 

6.5 

235 

3.175 

-23.9 

6.5 

236 

3.2 

-23.9 

6.5 

237 

3.225 

-23.9 

6.5 

238 

3.25 

-23.9 

6.5 

239 

0.1 

30.0 

6.6 

240 

3.225 

-23.9 

6.5 

241 

3.26 

-23.9 

6.5 

242 

3.265 

-23.9 

6.5 

243 

3.27 

-23.9 

6.5 

244 

3.275 

-23.9 

6.5 

245 

3.28 

-23.9 

6.5 

246 

3.285 

-23.9 

6.5 

247 

3.29 

-23.9 

6.5 

248 

3.295 

-23.9 

6.5 

249 

3.3 

-23.9 

6.5 

250 

3.5 

-23.9 

6.5 

251 

3.75 

-23.9 

6.5 

252 

4.0 

-23.9 

6.5 

253 

4.25 

-23.9 

6.5 

254 

4.5 

-23.9 

6.5 

255 

4.75 

-23.9 

6.5 

256 

5.0 

-23.9 

6.5 

257 

5.25 

-23.9 

6.5 

258 

4.775 

-23.9 

6.5 

259 

4.8 

-23.9 

6.5 

260 

4.825 

-23.9 

6.5 

261 

4.85 

-23.9 

6.5 

262 

4.875 

-23.9 

6.5 

263 

4.9 

-23.9 

6.5 

264 

4.925 

-23.9 

6.5 

265 

4.95 

-23.9 

6.5 

266 

4.975 

-23.9 

6.5 

267 

0.1 

0.0 

7.0 

268 

-2.15 

-30.38 

6.5 

269 

4.88 

-23.9 

6.5 

270 

4.885 

-23.9 

6.5 

271 

4.89 

-23.9 

6.5 

272 

4.895 

-23.9 

6.5 

273 

4.125 

-23.9 

6.5 
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Table  A.l  Summary  of  computer  runs  (finished) 


Case 

a(0) 

Pa 

la 

Vr 

274 

-2.775 

-31.44 

500.0 

6.5 

275 

-1.764 

-32.96 

500.0 

6.5 

276 

1.748 

-29.26 

500.0 

6.5 

277 

0.6263 

-24.0 

500.0 

6.5 

278 

-0.825 

-33.02 

500.0 

6.5 

279 

6.0 

-37.5 

500.0 

6.5 

6.0 

-34.602 

500.0 

6.5 

281 

6.0 

-33.451 

500.0 

6.5 

282 

6.0 

-32.661 

500.0 

6.5 

283 

6.0 

-32.023 

500.0 

6.5 

284 

6.0 

-31.466 

500.0 

6.5 

285 

6.0 

-30.956 

500.0 

6.5 

286 

6.0 

-30.472 

500.0 

6.5 

287 

6.0 

-29.528 

500.0 

6.5 

288 

6.0 

-29.044 

500.0 

6.5 

289 

6.0 

-28.534 

500.0 

6.5 

6.0 

-27.977 

500.0 

6.5 

291 

6.0 

-27.339 

500.0 

6.5 

292 

6.0 

-26.549 

500.0 

6.5 

293 

6.0 

-25.398 

500.0 

6.5 

294 

1.7255 

-26.549 

500.0 

6.5 

295 

2.30118 

-26.549 

500.0 

6.5 

296 

3.75 

-26.549 

500.0 

6.5 

297 

1.7255 

-25.398 

500.0 

6.5 

298 

2.30118 

-25.398 

500.0 

6.5 

299 

3.75 

-25.398 

500.0 

6.5 
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